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1 Introduction 



The evaluation of multi-loop Feynman diagrams has a long history dating back to the days of 
renormalization of QED (for an historical review see ref. [0]). In this respect there are theories 
more simple than others, noticeably QED and also to some extent QCD, where we have few masses 
and the analytical approach can be pushed very far. Typical examples are the calculations of g — 2 
in QED, see for instance ref. ||, or recent four-loop calculations in QCD 0]. Conversely the full 
electroweak Lagrangian shows several masses, ranging over a wide interval of values, therefore 
making the analytical evaluation of Feynman diagrams a complicated task. 

The modern era of loop-calculus begins with the work of 't Hooft and Veltman |J for arbitrary 
one- loop diagrams. As a result, an arbitrary scalar triangle diagram is expressed as the combina- 
tion of 12 di-logarithms (sometimes called Spence-functions) while 108 di-logarithms are needed 
for the general box diagram. 

In recent years we have been witnessing a huge amount of work in the direction of analytical 
evaluation of multi-loop diagrams. To be more precise we should distinguish between reduction 
procedures to bring specific classes of physics problems to some specific classes of integrals and 
calculational procedures for those specific classes of integrals. 

It would be impossible to quote the hundreds of papers that have been dealing with the 
subject and we are forced to make some selection. In the following we will briefly review the 
relevant literature, confining the search to those papers that are not directly connected to the 
present work. Other important references will appear in the body of the paper. 

Following and improving the techniques of || dimensionally regulated pentagon integrals have 
been considered in At the same time the methodology for one-loop n-point gauge theory 
amplitudes has been further developed in @ and in ||. 



The algebraic reduction of one loop graphs, introduced in ||, has been reconsidered in [|10 
More generally, the problem of reducing Feynman graph amplitudes to a minimal set of basic 
integrals has been discussed in fy] . 

An example of systematic approach to multi-loop calculations with the use of the program 
SCHOONSCHIP [0 is due to the Dubna group and a summary of the results can be found in 
ref. [jTB]. Meanwhile, algebraic programs like Mincer have been introduced |TJ 



An important breakthrough in this line of research is also represented by the algebraic approach, 
the so-called method of integration by parts. The simplest identity of this type was buried in a 



proof in [58]. A complete summary of everything relevant for 3-loop self-energies, as used in the 



next-to-next order (NNLO) calculations, is collected in [IS|. The concept and calculability of 
massless three-loop self-energies was reported in |TJJ. Further details for three- loop self-energies 
can be found in [|Tj]]. Also important in this context are explicit solutions of the recurrence 



relations, see for instance ref. 18 



Another line of development is related to the evaluation of massive integrals, e.g. in the context 
of effective heavy- mass theory. The most recent reference is |[L9| . Finally a group of people were 
successfully using the so-called method of uniqueness with various improvements, see ref. 0] and 
also JZT|,|21| and 
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Another significant development is represented by the method of differential equations, first 
discussed in []2l| and fully developed in |25[ (see also ref. p6[ ). Explicit solutions of the recurrence 
relations with respect to space-time dimension are shown in |27j . Generalized recurrence relations 
for two-loop propagator integrals with arbitrary masses are discussed in 



2s 



Furthermore, there is the subject of asymptotic expansions of Feynman diagrams. The first 
journal report of existence of simple formulas for OPE (those used in NNLO calculations) was in 
ref. ||29|| , while a summary of the rules was presented in ref. [[Kj. The euclidean variant of the 
theory appeared in ref. and ||39|| . The theory was reviewed - and the key trick to extend 
it to Minkowski space indicated - in ref. |32| . The final step was the extension of the theory 
to arbitrary diagrams in Minkowski space (including real phase-space diagrams), see for instance 



ref. 



Successive contributions can be found in ref. 
The techniques of Gegenbauer polynomials has proven to be useful for massless self-energies, 
mostly in two loops but with non-canonical powers of propagators, see ref. |37J and ref. B3] . 

Another meaningful technique, namely momentum expansion of multi-loop Feynman diagrams, 
was introduced in JIT]. Furthermore, the connection between Feynman integrals having different 
values of the space-time dimension was shown in |[42|, and the conformal mapping and Pade 



approximants to the calculation of various two-loop Feynman diagrams in 

Another significant attempt towards the analytical evaluation of multi-loop diagrams is repre- 
sented by the work of ref. [H - 

As far as applications to the standard model are concerned, two-loop self-energies have been 
discussed in [^0| . From the point of view of calculations with a direct relevance to standard model 
phenomenology we quote the next-to- leading two- loop electroweak top corrections of ref. [H5| , the 
fermionic two-loop contributions to muon decay and M w — M z interdependence of ref. [f|6| and 
the complete two- loop QED contributions to the muon lifetime |J47[| . 

All these approaches, with very few exceptions, have a common motivation: to compute an- 
alytically as much as possible of multi-loop Feynman diagrams. Soon or later this approach will 
collapse and one can easily foresee that this failure will show up at the level of a complete two- 
loop calculation in the standard model which is required to produce quantities as sin 2 9 l eS with a 
theoretical precision of 1 x 10~ 6 . For this reason one is lead to consider an alternative approach 
to the whole problem, namely to abandon the analytical way in favor of a fully automatic, nu- 
merical evaluation of multi-loop diagrams. The Holy Graal (General Recursive Applicative and 
Algorithmic Language) requires fast and accurate procedures to deal with the singularities of an 
arbitrary diagram. 

The outline of the paper will be as follows: in Section Q we briefly review the Bernstein-Tkachov 
(BT) theorem and in Section |2.1| we consider its application to one-loop diagrams. In particular, 
the infrared divergent triangle graph is discussed in Section p.l.lj , while other one-loop functions 
are studied in Section |2.1.2| and |2.1.3| . Reduction of tensor one-loop integrals is examined in Sec- 
tion |2.1.4j . Singularities of one-loop diagrams and the connection between generalized Bernstein 
functional relation and Landau singularities is briefly addressed in Section |2.1.5| . Multi-loop dia- 
grams are introduced in Section ^]2| and the main result of this paper, the minimal BT approach, 
is presented in Section 0. In Section O] and |3.1.3| we discuss the evaluation of the so-called sunset 
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two- loop diagram (or sunrise, depending on the mood of the author). Integration with complex 
Feynman parameters is introduced in Section |3.1.4| and a solution for numerical evaluation around 
and at the normal threshold is given in Section |3.1.5| . Special cases of the sunset diagram are 
discussed in Section |3.1.6| where numerical results are shown and comparison are performed with 
those results that are available in the literature. Tensor integrals are introduced and discussed in 
Section CO. Technical details are presented in several appendices. 



2 Bernstein-Tkachov theorem 

For a fast and accurate numerical evaluation of multi-loop diagrams there is some interesting 
proposal in the literature that has not yet received the due attention. 

The Bernstein-Tkachov theorem tells us that for any finite set of polynomials Vi(x), where 
vector of Feynman parameters, there exists an identity of the following form: 

V (x, d) J] Vt l+ \x) = 5 II Vi H (x). (1) 

i i 

where V is a polynomial of x and di = d/dxf, B and all coefficients of V are polynomials of fa 
and of the coefficients of Vi(x). The proof of the theorem uses methods of abstract algebra and 
we refer to the paper by F. V. Thachov |5(J for details. Given any Feynman diagram the fa in 



Eq.flip will be of the form —rii — s/2 where rtj are positive integers and n = 4 — e with n being the 
space-time dimension. One can apply recursively Eq.(|]) till the moment when all powers are of 
the form Ni—e/2 with N( > 0. The Laurent expansion in e yields the final form of the integrand. 



As pointed out in ref.[[5T|, with iVj = the imaginary part of the integrand is discontinuous, with 
Ni = 1 the integrand is continuous, with iVj = 2 has continuous first derivative and so on. In the 
next section we will show the solution for arbitrary one-loop diagrams. 

2.1 The one-loop case 

For arbitrary one-loop diagrams we have a universal master formula, again due to F. V. Tka- 
chov |5(J (see also ref. |51| and ref. f52|). Any one-loop Feynman diagram G, irrespectively from 



the number of external legs, is expressible as 

G= f dxV-»(x), (2) 

where the integration region is Xi > 0, J2i x % ^ 1 an d where V(x) is a quadratic form of x, 

V(x) =x l Hx + 2K l x + L. (3) 
The solution to the problem of determining the polynomial V is as follows: 

V = l- ( X + X ^ d ^ X = K t H-\ B = L-K t H- 1 K. (4) 
2(/i+l) 
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Figure 1: The three-point Green function. All external momenta are flowing inwards. 



Let us consider an example, the scalar 3-point function Cq of Fig. [T] The corresponding expression 
is written in terms of a 2-dimensional integral. We start from 



Co = C (pI,pI, (pi +P2) 2 ;mi,m 2 ,m 3 ) 



JL. 

i 7T 2 



d n q 



(q 2 + m\) ( (q + pif + m|J ( (<? + pi + P2) 2 + ' 



(5) 



where n = 4 — e is the space-time dimension and is the arbitrary unit of mass. After introducing 
Feynman parameters, we obtain 

C = dx dyV^ 1 ^ £ ^ 2 (x,y), V(x 1 y) = ax 2 + by 2 + cxy + dx + ey + f — id (6) 
Jo Jo 

where 5^0+ and where the coefficients of the quadratic form V are 



a = -Pi, = -Pi, c = pf + P2 - Q z , d = p A 2 + rri2- m'g, 

2 1 ,o2 1 2 2 r 2 

e = -P 2 + Q + mi - m 2 , / = m 3 . 



(7) 



Therefore, the quadratic form in Feynman parameters is specified by 



„ 1 (2a c 
H= 2 [c 2b 



K l = - (d, e) 



f-i5. 



The main idea beyond the BT-theorem is that we can now integrate by parts and raise powers in 
the integrand. The first step will be as follows: define (X, Y) = (H^^K and derive 



dx I dy V~ € ' 2 (x, y) + — / dx 
Jo Be Jo 



1 + Y) V~ e/2 (x,l) 
- Y V- e/2 (x, 0) + (1 + X) V^ t2 (1, x) - (x + X) V- e/2 



[X, X 



(9) 



Although this formulas could be used directly as it stands, better numerical convergency is reached 
when we raise again the power of V. This requires applying again Eq.(§) to the first term in Eq.(^) 
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with /i = — e. For the terms inside the second integral of Eq.@ we have now several quadratic 
forms in one variable to which the algorithm has to be applied: 



1) V(x, 1) = ax 2 + (c + d)x + b + e + f - id 

2) V(x, 0) = ax 2 + dx + / - iS 



3) V{l,x) 

4) V(x,x) 

The corresponding coefficients are 



1) H = 


a, 


K = 


2) H = 


a. 


K = 


3) H = 


b, 


K = 


4) H = 


a + b + c, 


K = 



bx 2 + (c + e)x + a + d + / - i5 
(a + b + c)x 2 + (d + e)x + f — id. 



Uc + d), L = b + e + f-i5; 
\d, L = f-i5; 



(10) 



i 



c + e), L = a + d + f — i5; 



The full result can be represented as the sum of a two-dimensional integral, a one-dimensional 
integral and a remainder, 



Gn = A 



-1 



A C* + G S + G U 



'12) 



Several combinations of the external parameters, not to be confused with Gram's determinants, 
enter into the final expression. They are 



G 2 = Aab — c 2 , Gii = a + b + c, G i2 
From them additional auxiliary quantities are constructed, 



a, G 



13 



b. 



(13) 



A = fG 2 - (bd 2 - cde + ae 2 ), 

Al = /GlI _(!*+£>!, a 2 = 



ci 2 



A 3 = (a + d + f), 



(14) 



and also 



a x = 2db — ce, 

d + e 
a>ix — — - — > 



= 2ea — dc, 
d c + e 

2> a 3x = 



(15) 



The three contributions of Eq . (|T2|) will then be written in terms of the auxiliary function V L (x, y) 
V(x, y) In V(x, y). We obtain 



dx 



dyV L (x,y), 



(16) 



Co = \ I 1 dx{v L (x, x)(a x - ay) [C^A" 1 + ^ G u A^ 1 ] + V L (x, 0)a y ^A" 1 + ^ G^A^ 1 

2 Jo 1 I J L I 
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/ 3 1 \ 3 

V L (l,x) [G 2 (A- 1 a x + - 13 Az 1 + -G 2 ) + -a x G 13 A 3 1 ]}, 



C ° = -~ G 2 A^ [a 3x (b + c + e) + i G 13 (^6 + c + e) 



12 2 v 



3 1 

b + -c + Ad + 2e + 6/) + -a^Ai (a + 6 + c + c/ + e) 



1 1 3 3 

- -a x a 3x A^ (b + c + e) + - a x G n A^ 1 (a + b + c + -d + -e) 
4 6 4 4 

14 1 

- - a x G 13 A^(-b + c + e) - - a 2/ a 1:c A^ 1 (a + 6 + c + c/ + e) 

o 4 

1 1 3 3 1 3 
+ ^ a y a 2x A2 1 (a + d) — - a y G u A^ 1 (a + b + c + -d + -e) + -aj / Gi 2 A 2 : (a + -d) 

+ - Vt(0, 0) [a^a^A^ 1 - a^a^A^ 1 + a y a 2x A2 1 



1 

4 
1 



Vt (1, 0) G 2 a3x A 3 1 + a x a 3x A 3 1 + 4 a y a 2x A 2 1 + a y G 12 A 2 1 



+ - 1^(1, 1) [G 2 a 3 xA 3 1 + G 2 Gi 3 A 3 1 - a x a lx A 1 1 + a^ag^Ag 1 



- a x G n A 1 1 + a x G 13 A 3 1 + a y a lx A 1 1 + a y G n A 1 



-1 



-1 



-1 



(17) 



18) 



It should be mentioned that in this approach the number of terms of the integrand really does 
not matter. The only relevant request to be made is that the algorithm is capable of produc- 
ing a smooth integrand. Furthermore the algorithm should be flexible enough to deal with all 
singularities, including infrared ones. 

2.1.1 Infrared divergent Co-function 

As a first example that this method can handle infrared divergences (IR) in the context of dimen- 
sional regularization we consider the case of the IR-divergent Co-function. From Eq.(f|) one can 
prove that B = when 



f = m 2 



s — 2 m , d = —2 m 



(19) 



i.e. for the case where an infrared divergency is developed. Here we introduced s — — (jp\ +p 2 ) 2 . 
Since the zeros of B correspond to leading Landau singularities (see Section 2.1.5 ) we recover the 
connection between infrared singularities and the more general class of Landau ones: as shown in 
Eq . fliDD , an infrared singularity for Cq manifests itself independently of s. As a consequence we 
have that the quadratic form V satisfies the following equation: 











1 + P 'K + P '* 



V» +1 (x,y) = 0, 



(20) 
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where the solution for the polynomial P is given by 

{1 - x\ —y} 
2(/i+l) 



P{x;y} 



(21) 



We use Eq.(|20|) to write 



dx / dyV- x - £te {x,y) 
o Jo 



f 1 dxV- x - ete { 
Jo 



e Jo 



(22) 



Since the IR pole at n = 4 will remain in the final answer we must keep trace of all n-dependent 
factors, 



Co (-m 2 , -m 2 , -s; m, 0, m) = - //7T e/2 Y (\ + - J dx V~ x ~ ele {x, 

= --fi £ - f 1 dxV- l - £/e {x,x), 
2 e Jo 



x) 



where we have introduced the usual ultraviolet and infrared regulators |>3 

12 . 12 



— 7 — hi7r, 



- + 7 + In 71, 



EE EE' 

and recall that these regulators satisfy the following relevant identity |53| : 

i + i = o. 

£ £ 



(23) 



(24) 



(25) 



Therefore, from Eq.(|2"3"D we recover the familiar result for the infrared Co-function, see ref []53 and 
also ref. ||. The power in V(x, x) can be raised once more with the help of the relation 



1 + 



2(^ + 1) V2 



- — x ) d. 



V 1+f *(x,x) = -- (s-4m 2 ) V^(x,x). 



(26) 



Collecting the various terms we obtain the final expression for the Co-function, 



G 



IR 



1 {-11 

s — 4 m 2 ^ 2 i 



If 1 ,, V(x,x) 
1 + - / dx\n I 



2 Jo 



1 
2 



l r 1 



1 , m 2 ( 1 , m 2N 
+ - In — 1 + - In — 

4 /x 2 V 2 /j, 2 I A Jo 



dx In 



^(a;, x) r 1 V^(x, x) 



3 + - In 



}. (27) 



In conclusion the BT-algorithm automatically extracts the poles in n — 4, both ultraviolet and 
infrared. 
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2.1.2 Other one-loop functions 



The procedure outlined above can be applied to any one-loop scalar function. Here we only show 
how it works for the massless pentagon (the massive case is as easy as the massless one). One can 
prove that 



E 



l—Xl 



1—XX—X2 



dx? 



1—XX—X2-X3 



dx& 



x l H~ 1 x + 2K t x 



-3-e/2 



(2f 



where the 4x4 matrix H is defined in terms of Mandelstam invariants, s^- = — (pi +Pj) 2 , by 
the following form: 



H 



1 
2 



/ — S51 S12 — S34 S45 ^ 

2 s 5 i -S34 - s 5 i s 23 - s 5 i 

- - -2s 34 -S34 

V- - - J 



and the vector K becomes 

K x = 0, K 2 



S51, #3 



S34, 



^4 = 0, 



(29) 



Starting from Eq.(^) we raise the power from —3 — e/2 to 1 — e/2 and, successively, integrate by 
parts. Note that, for the massless pentagon the coefficient B becomes, 



B — — S12S23S34S45S51. 
Id 



(30) 



For the fully massive pentagon we derive, after four iterations of the BT-algorithm (from power 
— 3— e/2 to power 1— e/2), an expression which, remarkably enough, has no 4-dimensional integrals. 
In Tab. [I] we give the number of terms for the scalar box, pentagon and hexagon functions, where 
in all cases we have raised powers up to 1 — e/2 an then we have expanded around e = up to 
terms of O (1). Once more, the number of terms is not the relevant issue, only the smoothness of 
the integrand matters. No attempt has been performed, so far, towards a systematic study of the 
singularities (zeros of B) of multi-leg, n > 5, one-loop diagrams. 



2.1.3 Other one-loop infrared-divergent functions 

One of the nice features of the method is that the coefficient B in Eq.([L]) contains all divergences 
of the diagram, including the infrared pole, as discussed for the triangle one-loop diagram in 
Section [2.1.1| . For other one-loop diagrams there is no need to have separate derivations for the 
infrared cases. All of them are reducible to Co-functions. Typical example is the box diagram Do, 
where we can use the following decomposition || (see also ref. 



A) (-ml, -m 2 e , -mj, -m 2 f , Q 2 , P 2 ; 0, m e , 0, m/) = ^ - J 77 (Gil 2 , P 2 ; m e , 



m 
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diagram 


5-dim 


4- dim 


3- dim 


2-dim 


1-dim 


0-dim 


D Box 






2 


54 


217 


98 


E Pentagon 






9 


189 


945 


12826 


F Hexagon 


1 


44 


693 


4581 


13860 


228846 



Table 1: Number of terms in scalar box, pentagon and hexagon functions where the BT-algorithm has 
been applied repeatedly up to powers 1 — e/2. 



+C {-ml, -m 2 f , P 2 ; m e , 0, m/j + C ~ m l, P 2 ', m f, 0, m e ) 



D {—ml, —ml, —m 



2 f ,-m 2 f ,Q 2 ,P 2 ;0,m e ,M z ,m f ) = ^ 2 ^ M2 [- J 7Z {Q 2 , P 2 ;m e ,m f 



-C {-m 2 e , -m 2 f ,P 2 ;m e ,M z ,m/) + C -m 2 e , P 2 ;m f ,0,m e ^ 

D {-m 2 e , -m 2 e , -m), -m),Q 2 , P 2 ; M z ,m e ,0,m f ^j = g 2 + M 2 [^t {Q 2 ,P 2 ;m e ,m f ^j 

+C {-m 2 e , -m 2 f ,P 2 ;m e ,0,m/) - C {~m 2 , -m 2 e ,P 2 ;m f , M z ,m^ . 

All infrared divergent terms are contained in the Co-functions and the additional integrals are as 
follows: 



2q-(q + Q) 



m 2 J 77 {Q 2 ,P 2 ;m e ,m f ) = /i 4 "" J d n q 

^ ( Q2 ' me ' mf ) = ^ I dnq d (0) dl (m e )d 2 (M z )d 3 (m f ) ' 
trc 2 J Z1 {Q 2 , P 2 ; m e , m f ) = /i 4 "" J d n q- 2Q ' {<1 + Q) 



d (0) di (m e ) d 2 (0) d 3 (to/) ' 
2q-Q 



(31) 



The propagators are defined by 
do = q 2 + m\ — i8, 



l do{M z )d l (m e )d 2 (0) d 3 (m f ) ' 



di — (q + Pi) 2 + to 2 , - i5, 



d 2 = (q + Pi + p 2 ) 2 + mj- iS, d 3 = (q + p 1 + p 2 + p 3 ) 2 + m\- id, 
with all four-momenta flowing inwards and we have explicitly written the internal masses. 



(32) 



9 



2.1.4 Reduction 



The standard procedure in dealing with one-loop diagrams is based on the introduction of scalar 
functions, -Bo, Co etc, and of tensor integrals that are reduced to linear combinations of scalar 
ones |J. The functional relation, Eq.(JI|), can be written in a more general form starting from 

/ dxQ{x)l[Vt(x), (33) 

where Q(x) is an arbitrary polynomial of x. Hence we are led to consider also the tensorial integrals 
in parametric space with no pre-reduction in momentum space and no appearance of standard 
Gram's determinants. Consider some simple example: we begin with the one-loop electron self- 
energy in QED, see Fig. ||, 




Figure 2: Fermionic self-energy in QED. 



x{x) = (1 — x) {xp 2 + m 2 ) 
After integration by parts we obtain 
S 



dx 



i (2 — e) xj> + (4 — e) m e 



X 



-e/2 





in 2 e 2 EpZ^-F 


1 

e 


p 2 

+ f 

(p 2 + m 2 e ) 


+ 1 


dx \ 

^2 2 


4 


p 2 


£ 


H 

(p 2 + m 2 ) 



2\ 2 



7 (p 2 + m 2 e ) 
2 p 2 

- — Ax) x In 



1 2 1 2 

- m„ p 

3 e Q 1 



/i 2 



27 -In 



1\ (p 2 + m 2 e f_ + p2 ln m 2 e 



2 ~3 P 



+ m z e ^3 In -f 



6 



dx x In ^7 
/i 2 



(34) 



(35) 



Note that the general idea is now different from the one which is at the basis of the analytical 
approach. The number of terms in the final expression does not matter and the only relevant 
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requisite is the smoothness of the integrand accompanied by the absence of high negative powers of 
true Gram's determinants in the reduction procedure. However, in this case a new and unpleasant 
feature shows up, namely the evaluation of the self-energy is numerical unstable at threshold, 

2 2 

p « — m e . 

Another example is given by the one-loop vertex of Fig. || which is reducible to the following 
expression: 



A, 



(2n) % 



For V\ we obtain 



16tt 2 



,2„3 



7/i 



V 1 (g 2 ;m,m) + {pi+p 2 ) u V 2 (Q 5 



(36) 



dx 



dy { [-2 ml V? - 2 Q 2 V?} x' 1 - ^ inS }, 



V7 = x 2 + 3 x - 3, V? = y 2 - xy + x - 1, 

vf ng 



I (e ~ 2) 2 



y) = x 2 m 2 e +y (x-y) Q l . 



(37) 



Note that we have already expanded in e because raising of powers can actually be done even in 
four dimensions. This new result can be derived by the following example: for /i = — 1 + e we 
write 



1 - 



(x + X) d x 



V e = BV 



-l+e 



(38) 



2 (ji+1) 

and expand in e before integration by parts. By equating the coefficients of e n we derive, for 
instance, that 



BV' 1 = 1- - {x + X) d x In V, etc, 



(39) 



showing that we can expand first and then integrate by parts. A similar result can be derived 
with n — e, giving 

B In V = V In V + ^ (x + X) d x [V (1 - In V)] . (40) 
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2.2.5 Zeros of B 



The relation between the zeros of B and the Landau singularities for an arbitrary diagram G 
remains to be examined, although B = should contain all singularities of G. In this section we 
consider the one-loop three-point function and establish the result that the coefficient B emerging 
after the first application of Eq.(JJ) is connected with the leading Landau singularity, the so- 
called anomalous threshold. Moreover, a second iteration produces several coefficients that are 
connected with sub-leading singularities, i.e. those of the contracted diagrams. This result holds 
also for higher one-loop functions, although the determination of the explicit form of the anomalous 
thresholds is, in general, an arduous task. 

Note that B is not a Gram determinant but its zeros are a vexation, as much as those of 
Gram determinants in the standard reduction procedure ||. As already observed in ref. the 
numerical advantage of smoother integrands is lost near thresholds where numerical convergence 
may become poor. This fact is best illustrated with a simple example, the scalar two point 
function, 

£ (- s;m)m ) = I-ln^ + 2-/31n^, (3 2 = I - A — . (41) 
s n p + 1 s 

The explicit result shows the well-known normal threshold at s = 4 m 2 . If we apply Eq.(|l]) the 
following result is obtained: 

B (s; m, m) = vr^ 2 // T (0 jf * dx X ^ 2 = j 2 vr"^ 2 // ^ [(3 - e) 



rl 

X / (/.r \ 



(42) 



with x — sx2 — sx + m 2 —i8. Of course there is no pole at s = 4 m 2 but numerical convergence 
becomes questionable. An interesting question is related to the singularity of B . Can we find it 
without having to perform the integration? The integrand is ln(x 2 — x + fj, 2 ) with /z 2 = m 2 /s and 
shows two branch points, x±, that are complex conjugated below threshold and pinch the contour 
at x — 1/2 when we are at threshold. We rewrite 

f 1 dx lux = - 2 + ln^ 2 + - (a^i 2 - l) ^dxx' 1 - (43) 
Jo 2 v ' Jo 

The integrand will the show poles at x± and we are more familiar with the fact that two poles, 



pinching the integration contour, produce a singularity for the function [49]. In ref. pOf an 



alternative method is suggested to deal with regions near thresholds (see also ref. |p4| ), while 



our solution to the problem will be discussed in Section |3.1.4j . Another example is represented by 
the triangle diagram. For the function 

C (-m 2 , -ml, -s;m fe ,M,m a ) (44) 
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we obtain the explicit form of the A factors of Eq. fll4|) : 



A x = -~ A (s, m 2 a , m 2 b ) , A 2 = -~ A (ikf 2 , m 2 , m 2 ) , A 3 = ~ A (m 2 , m 2 , 



I • • \ •< / - — •' | " \ -ai m b) i (45) 

that correspond to the three possible cuts of the triangle diagram, A being the Kallen function, 

\(x, y, z) = x 2 + y 2 + z 2 — 2 (xy + xz + yz) . (46) 

The zeros of Aj correspond to non-leading Landau singularities, a) normal thresholds at s — 
(m a + rrib) 2 etc. and b) pseudo-thresholds at s — (m a — mb) 2 etc. The latter do not show up in 
the physical Riemann sheet. Furthermore, we also derive 



A 



mi 



-M 2 s 



s + M 2 - 2 (m 2 + m 



For the special case m a = = m we get 

A = -M 2 s (s + M 2 — 4 



in 



(47) 



(4J 



and s = 4 m 2 — M 2 corresponds to the leading Landau singularity, the anomalous threshold. Note 
that the general form of the leading Landau singularity is 

s = -^hm 2 2 {ml +m 2 ) - {m\ + m 2 - M 2 ) (m 2 + m 2 -M 2 ) ± |~A(M 2 , m 2 , m 2 )A(M 2 , m 2 , m 2 )l V2 }, 

where we have introduced p 2 = —M 2 . 

For all the zeros of A and Aj numerical convergence is at stake. Here the method suffers the 
same disease of the usual appearance of Gram's determinants in the standard evaluation of Co and 
in the reduction of vector and tensor triangles, although the origin of the phenomenon is rather 
different. Finally, consider the simple case where m a = m& = 0. The exact expression is 



C (0, 0, -s; 0, M, 0) 



«2>-Li 2 (l + ]i L) 



(50) 



where ( is the Riemann zeta-function and Li 2 (z) is the standard di-logarithm. The zero of Gn is 
at s — 0, while G 2 has two zeros, respectively at s = and s = —M 2 . Note that there is no pole 
at s = in Eq. (|50D , that the function is regular at s = —M 2 and that the only singularity is the 
branch point of the di-logarithm. 

Actually these problems are easily solved in the standard approach. Consider the case of an 
arbitrary triangular graph for vanishing Gram's determinant: 



C (p 2 1 ,p 2 2J Q 2 ;m 1 ,m 2 ,m 3 ) = C = ^ C 0n (p 2 lJ p 2 2 ,Q 2 ;m 1 ,m 2j m 3 ) 8%. 



(51) 



n=0 
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The kinematics for vanishing Gram's determinant, A3 — > 0, is as follows: 



A 3 = P y 2 -( Pl . P2 y = 5,p\ 

a 2 n + S 3 ) P 



P 2 



Q 2 



(P1+P2) 



(1 - « ) 2 + #3 



Pi -P2 

P 2 . 



-a P 2 



(52) 



In lowest order the solution is to introduce three different pinches of the basic three-point integral 
which correspond to the three ways of cutting the triangular graph: 



a 



(i) 
00 



Bo (23) 



s 3 =o, 



a 



(2) 
00 



S (13) 



Then a solution, free of singularities, is 



G 



00 



1 

iVo 



Cg} + (l-a ) C&'+aoQ 



<5 3 =0, 



(2) 



a 



(3) 



00 



S (12) 



<5 3 =0 



,(1) 
00 



N = m 2 — m 3 + (m 1 — m 2 j oiq + P «o (1 — cxo) 
where ckq is given in Eq. (|52"D . At next order we have 12 terms: 



01 



-V 2 + I AiV - 3 C 00 + iV - 2 ^oo - 'i-^N n ~ 2 An (1) 



l^N - 2 A (2) + hxulN^Av (3) 



3 1 — «o 
1 

~ 3 



"0 ^0 



bi 2 N 2 Cqo ) 



01 iV °00 + o ai27V °00 _ o 7V U 00 

3 1 — «o 3 3 «o 

-ir»( 2 ) 1 ™ /vr-i^C 1 ) 



+ (1 - a ) iV- 1 ^ + a #o 
where we have introduced the following additional quantities: 



-fx 2 P 2 , 612 = 1 + n\ - l4, <& = 4//* + A 2 , 



A 2 = b 2 12 -Afxl 



Pi = — 
a 



1 



1 - a 



(53) 



(54) 



(55) 



(56) 



At second order we have an expression with 48 terms while 135 terms are needed at third order. 
For a general treatment of the reduction we refer to the work in ref. ||55||. 



2.2 Multi-loop diagrams 

In moving to multi-loop diagrams we recall that any diagram G with Nl legs and n\ loops is 
represent able as |56] 



G=^ 2 ) TH Y{N L -ln) J 



dx r , 8 (1 — 2V 



f/n/2 (y_ i5 jN L -nm/2: 



(57) 
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where T is the Euler gamma-function and where the integration measure can be written as 



N t TV, 

dx G = ]J dx i> x g = 52 x %- (58) 
Furthermore, the polynomials V and U are defined by 

v = J2 m2 x i + J2 x i - tj J2 B n * • * 

i i ij 

U = ^ II X * = d<3t C^™) ) f^rs = 51 X iVirVis, (59) 
T ir t £T i 

where r/j S is the projection of line i along the loop s. Furthermore, T is a co-tree and are 
the parametric functions for the given diagram. Although these functions can be determined 
completely by the topological structure of the diagram G we give a practical example of how to 
construct U and V for the two-loop diagram 5*4 of Fig. |] 




V2+P 



Figure 4: The two- loop diagram S4, Arrows indicate the momentum flow. 

After introducing Feynman parameters the integrand for S4 contains a factor l/D 4 with 

4 

D = "}2 Xi (qj t + ml) , qi = n, g 2 = r x - r 2 , g 3 = r 2 , g 4 = r 2 + p, (60) 

i=l 

where r s is the independent integration momentum around the loop s, Xi are Feynman parameters 
with J2i Xi = 1. The part of D which is quadratic in ri )2 will be written as 

r*Z7r, [/n=xi+x 2 , C/ 22 = x 2 + x 3 + x 4 , U 12 = U 2 i = —x 2 . (61) 

Next we rewrite C/jj as a sum, 

4 

Uij ViiVij x h (62) 

z=i 

and derive the coefficients rj as 

?7ll = ^21 = ??32 = ??4 2 = +1, V22 = -1, V3i = V4:i = Vu = 0. (63) 
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Furthermore, let U be the determinant of the matrix U^, thus 

U = det (Uij) = xi x 2 34 + x 2 £34, 



(64) 



where = Xi + Xj + . . . + x\. Momenta Pi will then be defined with p 4 = p and Pi = for i < 4. 
The following change of variables in the ri-integral is then performed: 



4 2 



1 



j=l t=l t=l 

Similarly we change variable also in the r2-integral, 

4 2 



<"EE 2^4^ ([/ 1 ) u = r^-^2 x 4 p^7] 4t (U = rf - x 4 jjrp 



u 



E E x jPjVjr (u 1 

j=l r=l 



^12 



We derive the following result: 



E x i (q 2 i + mf) ^r l Ur + V. 



which defines the polynomial V as 



(65) 



(66) 



(67) 



t r 2 1 2 L 2 2 

= 2^ x i m i +%4P -JjX l2 X 4 p. 

i=l ^ 

After a diagonalization of the symmetric matrix U, 

E(^ 1 )u> Ui'ji Aj/j = U Sij, 



(6£ 



(69) 



we perform a change of variables with unit Jacobian, Si = J2j ^ij r jy an d use 



n ^ 2 e^? + y l = / n + e + v 

i=l i=l ^ i=2 i=2 

W2-JV L 



-N L 



/2 n/ 2 T(N L -n/2) 
** Ul T(N L ) 



n *»< [E u ^ + v 



etc., 



i=2 i=2 

to obtain the result of Eq . (|57|) . Similarly for the diagram S3 of Fig. |7| we derive 

4 1 

U = Xi x 23 + x 2 x 3 , V = E m i + 77 ^1^2X3 P 2 - 



(70) 



(71) 



Note that UV-singularities come from U so that, for finite diagrams, one should raise only the 
factor 

V = U [E "i 2 Xi + E % 2 a?i] - E % ' Qj x ^ x r ( 72 ) 
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There are several comments to be made before we can start trying to apply the BT-theorem to 
any multi-loop diagram. For two-loop diagrams U is a quadratic form in the Feynman parameters 
x and V is a cubic, so that we have to construct a BT functional relation for a quintic or higher 
polynomial. For Nl = 3, 4, we have that the corresponding diagram is defined to consist of V to a 
positive power and, in principle, it is enough to regularize G with some Ks operation (S G G) |56| . 
For N L > 5, we have that the diagrams contains V to a negative power, U to a positive power. 
Therefore, there is no UV- divergency and we could raise the power for the cubic. 

For self-energies at zero external momentum, p 2 = 0, the procedure is quite straightforward. 
Consider, for example, the case n\ = 2, Nl = 3. First of all, we perform a projective transforma- 
tion 



x i — AiUij ' A? u v ^ ~ — 2> (^3) 

to obtain V — U, where U is a quadratic form. The important point being that we know how to 
raise powers for any quadratic form. 

Finally, Nl — 2 (ni + 1) is a nice example of a single polynomial of degree n\ + 1. Of course, 
vacuum diagrams are also easy. At two-loop level we have only the diagram of Fig. |5] which 
becomes 




Figure 5: The two- loop vacuum diagram with arbitrary masses. 



G vac 



X 



(4tt) 
1 - 



E-4 



// r (e - 1) 



K l H' l K 
1 



[x 



dx 



dy 

U £/2 (x,y), 



(74) 



where U can easily be derived according to the general rules presented above. The same is obviously 
true for tadpoles. All of this, however, still does not represent a real multi-loop application of the 
BT-theorem. 

Additional problems are represented by the fact that V is an highly incomplete polynomial 
with few non-zero coefficients. V can be found via a direct study of linear systems of n equations 
(usually, n ^> 100) in m variables and very often the matrix of coefficients has rank < n, so the 
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system is generally impossible and we need to go to higher values of n, m, i.e. beyond n min = 
min{n < m}. A simple counting shows the following situation: given a compact representation 
for V, 

V = P n + P^A + P^didj + ... (75) 

four variables and a cubic V require n eq = 126(330, 715) and n var = 155(415, 871) for n = 2(1, 0) 
and first (second, third) derivatives. Simple examples already show that, for realistic polynomials 
V, one has to go beyond second order in derivatives. 



3 A minimal BT approach 

The original proposal by F. V. Thachov requires knowledge of the polynomial V of Eq.(TT]), 
iterative application of the functional relation Eq.(|I]) followed by integration by parts. After a 
Laurent expansion in e = 4 — n one can achieve an arbitrary degree of smoothness of the integrand. 
In principle, the denominators B will contain thousands of terms and may lead to large numerical 
instabilities around thresholds. 

We have not been able to derive any compact form for the polynomial V of Eq. ([[]), even for 
the simplest two-loop topology. Instead, we have adopted a different strategy aimed to deal with 
arbitrary two-loop diagrams. It represents a compromise based on the simple observation that we 
know how to apply the BT-iterative procedure for arbitrary one-loop diagrams. Therefore, given 
any two-loop diagram G (for an illustration see Fig. |6|) we apply the BT functional relation to 
G L , the one- loop sub-diagram of G which has the largest number of internal lines. In this way 
the integrand for G in x-space can be made smooth, a part from the factor B of Eq. ([!]) which is 
now a polynomial in x s , the Feynman parameters needed for the complementary sub-diagram of 
G with the smallest number of internal lines, G s . The sub-diagram G s , after integration over its 
momentum, becomes an additional - x 5 -dependent and with non-canonical power - propagator 
forGj,. 

This procedure can be immediately generalized to any number of loops. Furthermore, one 
should realize that the BT procedure does not introduce singularities through B(x s ), a part from 
the singularities in the external parameters of the original diagram. Therefore, before performing 
the x s -integration we move the integration contour into the complex hyper-plane, thus avoiding 
the crossing of apparent singularities. The idea of distorting the contour has been introduced 
in |5P[ and we have modified it to deal with the ^-integration. 

A complete study of two- loop two-point functions (self-energies) has shown |57j that, for one 



iteration, a transformation of the corresponding Feynman parameters can always be found that 
produces a coefficient B which is x s -independent. In this way we are able to avoid distortion 
of the integration contour. However, this coefficient B will vanish at some non-leading Landau 
singularity of the diagram where additional analytical work is needed before starting the numerical 
evaluation. A detailed derivation will be given in f57j. 

In the following section we present a simple example that illustrates all the features connected 
with integration in complex parametric space. 
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Figure 6: A generic two-loop diagram G; internal solid lines give the largest sub-diagram G L . 



3.0.1 An example 

An instructive example is the following: suppose that we have to compute 



/ dx B (— s; xm,m) , 
J o 



(76) 



where Bq is the scalar, one-loop, two-point function. Suppose also that we have no analytical 
expression available for Bq which, of course, we have but, nevertheless, let us assume that we do 
not. The integral of Eq. (|76|) becomes 



27 Jo 



dx / dyx 



-e/2 



I = i n 2 ~^ 2 T 
x(x, y) = sy 2 — (s — m 2 + x 2 m 2 ^ + x 2 m 2 — id, 
with 5 — > + . Then Eq.(H|) applies for 



As V 

with the following result: 



2 2 2 

. . ,s,x m ,m 
As • 



v / 2,22 

Y = [s — m + x m 

2s 



dy X ~ t/2 = -4 



After integration by parts we obtain 

s 1 



A Jo 



dy 



X 



l-e/2 



dy x ' e ' 2 = -A ~ ^— [(3 - e) jT dy x 1 ^ 2 - (1 + Y) (m 2 ) 1 ^ + Y (x 2 m 2 ) 



l-e/2- 



(77) 



(7S 



(79) 



*0) 



This strategy is dubious if the roots of A = are inside the interval [0, 1], on the real axis. These 
roots are easily evaluated and are given by 



x\ 



1±£ 



2 1 

, jj, s = m . 



Therefore 



x = — , for u > 1, 

fi 



x = , for n < -, 

H 2 



(82) 



19 



are the apparent singularities in the interval [0, 1]. Crossing of these points is avoided by moving 
the x-integration into the complex plane. How to deform the integration path is suggested by 
simple considerations: take the quadratic form \ of Eq.([77]); as a function of x there are two roots, 



X{x,y) = ^ 2 (1-2/) (x + xo + iS) (x-x -iS), x 2 = — ■ r — . (83) 

Therefore, for 1 — yU 2 < y < 1 we have two branch points at ±z|xq|- Instead, for < y < 1 — fj 2 
the branch points are real and x < f° r < x < Xo- Furthermore, Xq < 1 for y < y_ or y > y + 
with y± = (1 ± (3)/ 2 and f3 2 = 1 — 4/i 2 . Thanks to the infinitesimal imaginary parts we have 

In x = In fJ 2 + In (1 — y) + In (x + xq + i S) + In (x — xo — i 5) , (84) 

and only the last logarithm contributes to the imaginary part. To avoid crossing the cut of the 
logarithm ([— oo,xo]) we choose any smooth curve connecting x = to x = 1 in the x lower half- 
plane. By comparing two expressions for / (as defined in Eq . (|76|) ) , the first computed from the 
known analytical result for the .Bo-function and a second one obtained through the double-integral 
originating from Eq.([30|) we find agreement in 8-digits, as shown in Tab. |2|. 



^[GeV] 




Eq-© 


1.1 


1.12256370(1) + i 0.042810523(5) 


1.12256370(1) + i 0.042810523(11) 


10 


-2.52793882(1) + i 3.09949189(0) 


-2.52793882(3) + i 3.09949189(3) 


100 


-7.20895671(4) + i 3.14117375(0) 


-7.20895671(7) +i 3.14117375(7) 



Table 2: / from Eq.© or Eq.© with m = 1. 



3.1 The S3 or sunset topology 

After some preliminar considerations we are ready to start with a realistic example. Consider 
the simplest two-loop, two-point topology with three internal lines, the so-called sunset diagram 
illustrated in Fig. 0. Literature relevant to the sunset (sunrise) topology is assembled, for conve- 
nience, in ref. [Q. The corresponding integral will be 

vr 4 S 3 = ^ J d n gi d n q 2 (q\ + m\) ' ( ( Ql - q 2 + p) 2 + m 2 .) ^ (q 2 + m 2 ) . (85) 

Here we introduce a special notation, (a, m\ \ 7, rri2,p \ /3, m^) , for 

/i 2e J d n qi d n q 2 (q 2 + m 2 y a ( ( 9l - q 2 + p) 2 + m 2 )" 7 (g^ + m 2 )"^, (86) 
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and (fix . . . ^ \ v\ . . . Vj \ a, mi | 7, m 2 , p \ {3, m 3 ) for 

^ J d n gi d n q 2 gf 1 . . . . . . q? (g? + m?) ~° ( (fc - q 2 + pf + mj) (q 2 2 + mf) ~ fi . (87) 




Figure 7: The S3 topology for two-loop self-energy diagrams. 



The best way of dealing with the sunset integral is to introduce the so-called (uf3^) partial 



integration |58| to obtain 



S 3 = J2 m ? S 3i + S - 



3p- 



i=l 



With a suitable change of variables the four functions of Eq.(^) can be cast into the following 
form: 



ir 4 S : 



31 



-^1 J d n qi d n q 2 (q\ + 



(qi - Q2 -pf + m 2 2 ) [Q$ + ml 



(e - 1) 1 (1, m 3 I 1, m 2 , -p | 2, m x ) 



7l 4 S: 



32 



^-j J d n qi d n q 2 (q\ + mj) 1 ( (q x - q 2 + pf + mf) ' (q 2 2 + mi 



(e - 1) (l,mi I l,m 3 ,p | 2,m 2 ) 



7l 4 S: 



33 



J d n Ql d n q 2 (ql + miy 1 ( ( Ql - q 2 + pf + mj) ^ (q 2 2 + 



in. 



(e - 1) 1 (l,TOi I l,m 3 ,p I 2,m 2 ) 



tt 4 5 



3p 



2c 



_7 / d n qid n q 2 q2-p (qi + mfj {{qi - q2 + pf + mfj (q% + m\ 
(e - (0 I p I l,mi I l,m 3 ,p | 2,m 2 ) , 



39) 



where, as usual, n = 4 — e. Their evaluation is discussed in Section 3.1.2 and in Section 3.1.3. 
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3.1.1 Landau equations for S; 



33 



Before starting the evaluation of the sunset diagram it is important to know as much as possible 
about its singularities as a function of p 2 and of the internal masses. The corresponding Landau 
equations are 



«i (ql + mf) = 0, a 2 ((qi~ q 2 + p) 2 + m\ 



0, 



"3 {ql + ml) = 0, 



(90) 



and also 

aiqi^ + a 2 (qi - q 2 + p)^ = 0, -a 2 (qi - q 2 + p)n + 0! 3 q 2fl = 0. (91) 

The leading Landau singularity occurs for ctj ^ 0,Vi. We multiply the two equations Eq. (|9T|) by 
qiii, q2fj, and p M respectively This gives an homogeneous system of six equations. If all ctj are 
different from zero, the singularity will occur for 

i 1 r o o 



ol 



-mi 



Ql 



6 i7 12 — ~m 3 , q\ ■ q 2 
The Landau equations become as follows: 



m — m, — mo 



s + 2p- (q 1 -q 2 ) 



(92) 



q2 ■ P + - (s — ml — ml + ml 



qi ■ p — <?2 • P — + m i — m 2 + m 3 
qi ■ P on + (s + qi ■ p - q 2 ■ p) a 2 = 0, 

1 / o o 



0! 2 = 

Qti + 



Qi ' V - ^ ( s + m i ~ m l ~ m l) a 2 = 0, 



-Q2 -P - 2^ s 



m 



a;2 + 



<?i • P - Q2 ■ P - ^ + m l ~ m 2 + m l) 



a 3 = 0, 



— <7i ■ p + - (s + m 2 — wij — m| 



«2 — ^3^3 

s ~ qi • P + <?2 • P) ol 2 + q 2 ■ p a 3 = 0. 



0. 



(93) 



There are two compatibility conditions for the first three equations that can be used to get q± ■ p 
and q 2 ■ p, 



<n-P= T 



mi s — mf + (mi + m 2 y 



Q2-P 



mi + m 2 

Inserting these solutions into Eq.(p3|) we obtain that 



Oil 



m 2 



m% 



mi + m2 



1 m 2 

— «2, «3 = o — 3 

?71i 2 ?7l2 

is a solution of the original system if and only if 

s = (mi + m 2 ± ms) 2 . 



s + mo — (mi + m 2 ) 



mi - m 2 \ a 2 , 



(94) 



(95) 



(96) 



These are the leading singularities (plus permutations of masses), referred to as normal and pseudo 
thresholds. Note that pseudo-threshold does not correspond to a singularity on the physical 
Riemann sheet. 
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3.1.2 Evaluation of S: 



3i 



We will now compute S 33 , as given in Eq.(^), and the remaining S^-terms will follow from S 33 
with the corresponding permutations: 



31 with mi <-> m,3, p — p, 

32 with m 2 <-> m 3 . 

First a Feynman parameter x is introduced, and the gi-integration is performed giving 



(97) 



S 33 = iir 



-2-e/2 _/f 



rr(l — x) (g — p) 2 + xm\ + (1 — x)m\ 



-6/2 



(9* 



where we have introduced the Euler T-function. Next we change integration variable, q' = q — p, 
and introduce dependent mass, 



S 33 can be cast into the following form: 

33 = I TT ' ; 









-" 2 J 


(0 


/ dx 


x(l — x) 




Jo 







(1 — x)m\ + xm\ 
x(l — x) 



d n q- 



(99) 



q 2 + ml 



(g 2 + m 2)lW2 ( {q+p f + 



mi 



(100) 



According to some standard procedure we introduce a new Feynman parameter y and write 



S 33 = in 



-2-e/2 f 1 



X 



e- 1 V2 

d n q (q 2 + - 2 



dx [ dy x(l — x) ^ (1 — yfl 2 y 



r 



j 



g +2yp-q + y [p + m 3 ) + (1 - y) m x 



-3-e/2 



Let C/ be defined by 



*/ = 3/(1 - y)p 2 + yrri3 + (1 - y)m 2 - i<J, 



(101) 
(102) 



with 5 — > 0+, then S33 becomes 



^33 = — 7T 



1 /•! 



- ^— T ( m ,, I' dx I' dy - x)]-^ (1 - y)^ 
e - 1 T (1 + e/2) 7o 7o L V 7J V ' 



x y 



r (1 + e) (yV + ml) IT 1 - + (2 - f ) T (e) [/ 



(103) 



Next we introduce the Mandelstam variable s, 



P 



(104) 



23 



and scale U accordingly: 

U = sx, 



x = -y (i - y) + i4v + A ( 1 -y)~ i 



(105) 



For the rest of this discussion we assume that s > and will postpone to Appendix B the discussion 
of the case s < 0. 

There is no need to evaluate numerically the poles at e = 0(n = 4). We perform a Laurent 
expansion of the integrand around e = 0, after which the double pole is easy to compute and gives 



4 + A 

e 2 



2 

UV1 



^uv = 7 + In 7r + In — . 

/i 2 



(106) 



In this result 7 = 0.577216 is the Euler's constant and we have included in the definition of double- 
pole a bunch of constants. Also the single pole at e = 0(n = 4) can be computed analytically. The 
x — y integration is greatly simplified by the use of the following relation: 



i4-y 2 



, y d 

1 

x oy 



(107) 



With the help of Eq.( |107|) we find the following result for the single-pole: 

-~A UV + (2 1n/4-3) (A vv -^\, 
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where, again, a bunch of constants have been included into the definition of single-pole. A detailed 
derivation is presented in Appendix A. Continuing our calculation we write 



S: 
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r (c/2) 



e-l r(l + e/2) V/A 

.2 , ,.2 N 



dx 







/2 







dy x(l -x) (1 - y) t/2 y 



(2--) r(c) x" 



x r(l + e) (-r + ^J V 
The quadratic form x is defined by 

X = ify 2 + 2ir 2/ + L 

with coefficients 



1 / 2 2 

-1-/^3 + ^ 



L = /4 



(109) 



(110) 
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Our strategy will be to single out the y-integral and to raise the power of x- The coefficients 
relevant for this case are: 



B = nl - i (1 - 14 + K 

F = -i(l 

9. V 



(4 + l4 



4 

/U = — 1 



(112) 
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where A is the usual Kallen function. Therefore we obtain 

4 



X 



-l-e 



x 



(113) 



No spurious singularity at e = is introduced in the r.h.s. of Eq.( |113| ) and to avoid apparent poles 
we use a new variant of the basic functional relation, Eq. ([[]): 



X 



-l-c 



4 {l-i (V + Y) |-lnx-6 [ln X -| (y + F)|-ln 2 X + 0(e 2 )]}. (114) 



9y 



Note that no new singularity in the variable x is introduced by the procedure. However, numerical 
instabilities could be introduced in the ^-integration, especially if the process of raising powers is 
continued, resulting in a final expression with higher and higher negative powers of A. For this 
reason we are led to study the zeros of the Kallen function, i.e. A (1, /i 2 , /i 2 ) = 0. Clearly, if s < 
there are no real solutions for x. Otherwise solutions are given by 



04) : 



Consider now /z 2 as a function of the variable x, 



[1 - x)^i\ + x^l 
x (1 — x) 



(115) 



(116) 



The minimum, for /z 2 , occurs at x± = ± /ii). Only x + lies between and 1, corresponding 

to 

C<Ln = ^+^- ( 117 ) 

There are three distinct possibilities: 

1. the root (/i 2 )+ is below the minimum, i.e. /ii + \ii — /13 > 1, therefore A can never be zero. 

2. Only one root is above the minimum, i.e. (1 — /13) 2 < (//i + /12) 2 < (1 + /is) 2 , when there 
are two values of x where A = 0, 



x 



L "2(1 + /i 3 ) 2 



;i + /^s) 2 + A - A ± a i/2 ( (i + /i 3 ) 2 , /i 2 , iA 
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3. Both roots are above the minimum, i.e. (/ii + /Z2) 2 < (1 — /is) 2 , when we have four values 
of x where A = 0. The new pair of points is given by 



x 4 



2 (i-/i 3 y 



^ - /^s) 2 + /i? - /i 2 ± A 1/2 ( (1 - /i 3 ) 2 , /i 2 , Z^) 



(119) 
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For the finite part of this Feynman diagram we have now to understand when and where the 
imaginary part develops. Let us assume that s > 0. Imaginary parts are present as soon as x < 0. 
There are two solutions for y to the equation \ = 0, 

y ± = ±[l-»l + fJ *± A 1/2 (l, fil /i 2 )] . (120) 

Therefore, lnx develops an imaginary part for y_ < y < y + , if and only if A(l, jti§, f/%.) > 0. As 
before, there are only three possibilities, 

1. (1 + ji\) 2 < (/ii + H2) 2 , where A is always positive, 

2. (1 — /i 3 ) 2 < (fii + fi 2 ) 2 < (1 + /U3) 2 , where we have two alternatives, 

(/ii + /i 2 ) 2 < /i 2 , < (1 + /i 3 ) 2 , A is negative, 

/-t 2 > (1 + , A is positive, (121) 

3. (/ii + yU 2 ) 2 < (1 — Z^ 2 ) 2 where we have three alternatives, 

(yUi + /U2) 2 < Ai 2 < (1 - /J-3) 2 , A is positive, 
1 — ) < ^ x < (1 + ^3) , A is negative, 

/-t 2 > (1 + /W3) 2 , A is positive, (122) 

Consider the original integral, Eq.( |109|) ; the integration hyper-contour can be distorted away 
from its original real location and when the possibility of this distortion ceases, we encounter a 
singularity of the function The difficulty in locating the singularities lies in imagining what 
happens in the multi-dimensional (complex) space of integration. In the case of Eq. (|109|) the 
integrand shows poles at the zeros of x{ x iU)i i- e - f° r f° r U = U±{ x ) with y± given in Eq.( |120| ). 
We may distort the interval y G [0,1] but a pinch may appear for y + = y_, which requires 
A(l, /i 2 , /i 2 ,) = 0. Again, we can move the x-integration contour to avoid the points where this 
happens, i.e. for /i 2 = (1 ± /i3) 2 , 

Hi = (1 ± /i 3 ) 2 , for x = x% (123) 

However, if x\ = x± or = xZ a singularity may appear. Note that for // 2 = (1 + /i 3 ) 2 we will 
have 

y± = i + Ms, x = (y-i-/i 3 ) 2 , (124) 

while for /i 2 = (1 — /i 3 ) 2 we have 

y± = 1-/^3, x = (2/- 1 + A i 3) 2 , (125) 

so that, in the first case, the pinch should correspond to y > 1, therefore there is no singularity. In 
this case we will speak of the so-called pseudo-threshold, corresponding to s — (mi + m 2 — m 3 ) 2 . 
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It is a well known fact that the S3 topology is regular at any of the pseudo-thresholds. On the 
contrary, when \j? x — > (1 — /i 3 ) 2 (or x — > x±) we have y + — > y_ from opposite sides of the integration 
contour but the relevant part of the integrand is (/^ — y 2 )/x which, for y = 1 — fi 3 + A (A — ► 0) 
and fj%, = (1 — /U 3 ) 2 behaves as 

ti-y 2 2 (1 - /i 3 ) a 

— t A*-i6 ■ (126) 

Therefore there is no pinch, even in this situation. Note that the condition x+ = xZ corresponds 
to the so-called normal threshold, s = (mi + 7712 + m 3 ) 2 which is, nevertheless, a singular point of 
the sunset diagram. The singularity clearly arises from the remaining x~ € terms of the integrand, 
in the same way as the normal threshold singularity arises in the one-loop self-energy. This rather 
long discussion has been introduced for one specific purpose: if we rewrite S33 as 

then ^33 (e, x) has no poles at x = x±. In other words points of analyticity for 

jf dy (1 - yf 2 y [v (1 + e) (-y 2 + ^) ^ + (2 - £) V (e) 
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This fact remains obviously true after raising powers of \i operation that brings negative powers 
of A in the result and, due to analyticity of the integrand in those points, we can distort the 
integration contour. One can reach this conclusion also by following an alternative derivation 
of the result. Starting from Eq. ( |109| ) we apply Eq. (|114j) only once. The integral will be of the 
following form: 

^ 3 = Jo dX Jo dV (X + J °) ' A = A(1 '^)- ( 129 ) 
A simple exercise shows that 

t dyh = 0, for /i 2 = (1 - ^f, (130) 
Jo 

so that we can replace 

h{x, y) ^ I suh (x, y) = h(x, y) ~h^ x = l- /i 3 ) . (131) 

Just one iteration of the raising procedure gives Jo,i continuous, as compared with two iterations 
that give also continuous first derivatives but, in this way, the whole integrand is well behaved 
around the zeros of A(l, /1 2 ). As for the one-loop, two-point function, once a suitable integral 
representation is found, the branch point corresponding to normal threshold is not due to poles 
of the integrand. 

We repeat that distortion of the contour is only dictated by the need of achieving numerical 
convergence and must satisfy one criterion: the new contour cannot cross the cuts of x~ e - Since 
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at the normal threshold the two branch points approach each other on the real axis, there we 
can no longer distort the contour and face numerical instabilities. Therefore, the method cannot 
be applied, as it stands, when we are exactly at threshold and an alternative algorithm will be 
presented in Section |3.1.5 . 

Summarizing, me may say that BT method has the advantage of moving poles into branch 
points and minimal BT method requires a knowledge of the analytical properties of a one-loop 
sub-diagram of the whole diagram G, a task much simpler than the one of knowing the whole 
analytical structure of G. 

After a discussion on the singularities of A and of the imaginary parts we continue by performing 
an integration by parts after the first operation of raising the powers. This will be achieved, after 
inserting Eq. (|114j) into Eq. (|109|) , with the help of the following relations: 



1 dyy n (l-y) e/2 |- \n m X (x,y) =-6 n0 ln X (x,0) 
dy 



dy 



ny 



n-l 



(1-1/) 



e/2 



\ia m x(x,y), 



[i-yf 2 - 1 

followed by a subtraction procedure that takes care of the 1 — y terms, 

f 1 dyy n (1 - yf 2 ~ l ln m X (x, y) = - ln m X (x, 1) + /' dy (1 - y)^ 1 
Jo e Jo 



(132) 



y n \n m X (x,y)-\n m x(x,l) 



The procedure of raising the power of \ can be repeated. Since 

X (x,y)=y 2 + 2Ky + L, 

we could use 



Bx 



1 



2 fl -el 



(y + Y) 



d_ 

dy 



X 



l-e 



(133) 



(134) 



(135) 



but it is better to expand this relation before integrating by parts, thus avoiding the introduction 
of fictitious poles at e = 0. The result that we obtain by expanding in powers of e and by 
subsequently equating coefficients of the same powers in e is as follows: 



d r 



B lux = X + ^ (y + Y) ^- [x (1 - lux) 

Q _ / 

B In 2 x = X In 2 X ~ {y + Y) ^- |x [1 - ln X H 



ln 2 x 



(136) 



In this way terms as ln x or In 2 x ar e transformed into x l n X an d X In 2 X- However, some care is 
needed in the presence of terms containing ln(l — y). In this case we proceed as follows: 

f 1 dyy n In (1-y) \n X = C dyy n ln (1 - y) L^ + \nA 

Jo Jo v f4 J 
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-^ln,i|;i + /^».n(l- !/ )ln^ 



(137) 



where use has been made of the relation x(x, 1) = /if. In this case the quadratic form in y to 
which we apply the algorithm is 



-2+2-2 y + -o, 



giving B and Y coefficients, 



pi 



^3 ^3 ^3 

Y = K. 

Therefore we will have terms of the following form: 

f 1 d 

/ dyy n In (1 - y) — (x lux) 

Jo ay 

which, after integration by parts, become 



1 x r — 1 y n 

dyxln— ny n ln(l-y)H - 

o f4 1 y-l 



:i38) 



(139) 



(140) 



(141) 



In this way we arrive at our final formulas for 5*33 which is shown explicitly in Appendix C. The 
other S 3i terms follow through the permutations of Eq. fl9"T|) . 



3.1.3 Evaluation of S; 



3p 



The remaining contribution is Ss p , defined in Eq.fl59p. Everything proceeds as in Section |3.1.2j , in 
particular the double-pole turns out to be zero and the single-pole is 



(142) 



with A uv defined in Eq. (|106|) . The complete expression for Ss p is also given in Appendix C. 
Collecting the various terms we obtain the following expression for 5*3: 



S, 



3 double- 



-pole = E m \ (i + A2 uv) 
i=l ve 7 



A r , v - 



'3 single— pole 



o \ "* 2 ^UV 



i=l 



+ ( A --~) [X> 2 ( 21n v- 3 ) + S- (143) 



The finite part will also be given in Appendix C. 
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3.1.4 Integration in the complex plane 



The integrand for S3 is now a smooth function but for inverse powers of A whose zeros are known 
and their crossing should be avoided by moving the ^-integration from the interval [0, 1] to the 
complex plane. Technically speaking, when we are exactly at the normal threshold, our method 
is no better than any other numerical method. Away from it, the actual shape of the integration 
path follows (as discussed in Section |3.0.1| ) from a careful examination of the imaginary part of the 



logarithms appearing in the integrand. The following recipe will produce the correct imaginary 
part for 5*3. If (1 — < (fix + /i 2 ) 2 < (1 + /is) 2 we compute x±, while for (1 — /i 3 ) 2 > («i + /i 2 ) 2 
we compute both 2^ and x±. Next we will have to investigate the branch points of the logarithms. 
The argument of the logarithms is a quadratic form in x, 

y) = -y(y-\ + M 2 ) x 2 + [y (y - 1 + /i 2 ) + (1 - y) (/i 2 - /i 2 )] X + (1 - y) /i 2 . (144) 

If the corresponding discriminant is positive we will have two real roots, x LtR , otherwise we have 
a pair of complex conjugated roots. In the latter case we select a contour which starts at x = 
and bypasses x = x± or x = x±, x± in the upper half-plane to return to x — 1. The distorted 
contour must avoid any crossing of the cuts of the logarithm which will be, in this case, parallel 
to the imaginary x-axis. 

If the roots are real me must distinguish between two cases. Let us denote the quadratic in x 
of Eq. ( |144| ) as ax 2 + bx + c, with b 2 > 4 ac. If a is positive then the cut is between x L and x R , while 
for negative a the cut is [— 00, x L ] U [x R , +00]. Furthermore, if x = a + i ft, the imaginary part 
of the logarithm is ft(2aa + b). This consideration immediately tells us the sign of the imaginary 
part when x approaches the real axis on the cut. This sign is crucial when we distort the contour 
in the complex plane since, for x real and on the cut, the i5 prescription, 

ln(ax 2 + bx + c) — > ln(ax 2 + bx + c — i 5) , (145) 

gives —m for the imaginary part. To give an example, when the branch points are real and a is 
positive, any contour should start at x — 0, bypass x±, return to the real axis from above and for 
Hex < —bj (2a), leave the real axis in the lower half-plane for Rex > — 6/ (2a) and bypass x\ to 
return to x — 1. Typical examples are shown in Figs. ITOl. 



We now examine the behavior of the integrand close to the normal threshold. Consider Fig. |TT 
which corresponds to s = (mi + ni2 + m^) 2 for some choice of masses: for y below some value yo 
we have that x L , R are complex conjugated and their real part is internal to the interval [0, 1]. Just 
below yo the two imaginary parts go to infinity and above yo the two branch points become real 
but are external to [0, 1]. Furthermore, for some value y t h the two complex roots x L<R pinch the 
integration contour at x = xZ = xT. Therefore, any distortion of it will inevitably cross the cut, 
giving the wrong answer. In a word, we cannot distort the contour any longer. Note that there is 
no crossing of poles implied here, only a serious problem of numerical convergence of the minimal 
BT approach due to negative powers of A(l, /x 2 ) which exhibits a double zero at x — xZ = 
For different choices of the internal masses, at threshold, we may have a different behavior of x L R 
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but there will always be a y t h where they are complex conjugated, with their real part approaching 



x 



x 



and with imaginary parts pinching the real axis. 



The situation at pseudo-threshold, s 



(mi + m 2 - m 3 j 



is different as one can see from 



Fig. [12]. Below some value yo the roots x LtR are complex conjugated but the absolute value of the 
imaginary part has a minimum and, therefore no pinch will occur. Above yo the roots are real but 
external to the interval [0, 1]. Therefore a distortion of the integration contour, to avoid the point 



x 



is always possible. 



3.1.5 The normal threshold for the sunset 



The arguments developed in Section j.l.4| show that we have to use some other algorithm at 
the normal threshold s = (mi + m 2 + m 3 ) 2 . Actually, for and S^ p we have a much simpler 
derivation of the result which, although different from the BT-approach, allows us to raise powers. 
We consider again Eq.( |109|) and apply directly Eq. ( |107j ). After the expansion around n = 4 we 
have the following integrations by parts: 



o 



dyy n ln(l -y)8 v \nx(x,y) 



1 , , x{x,y) 

dy In— — — 



ny 



n-l 



ln(l - y) 



y 



l-y J 

,n-l 



1 1 r f 1 

dy y n In x(x, y) d y In x(x, y) = - In 2 x(x, 1) - 5 n>0 In 2 x(x,0)-n / dy y n ~ x In 2 x(x, y) 
A 1 Jo J 

1 j-l 

dyy n d y \nx(x,y) =\nx(x,l) -5 nfi \nx(x,0) -n dyy 71 ' 1 In x(x,y) . (146) 
o Jo J 

Using these relations we obtain a very simple expression for the finite part of 5*33, 



S%» = jT 1 dx jf 1 dy [lnt(x,y) + 

r 1 r 1 3 
S'gp = dx dy,y\n£(x,y) + -, 
Jo Jo 8 



+ 1U/. 2 (ln/, 2 -4) + ^ + ic(2), 



(147) 



where £ = x(l — x)x and the '-(-'-distribution is defined by Eq.( |181| ) of Appendix C. Therefore, we 
have reached our goal of expressing the result in terms of a smooth integrand, without obnoxious 
overall factors that can vanish in certain regions with a consequent numerical instability. The 
reason why we have not indicated Eq.( 147 ) as the main result of our work is based on the fact 
that the same result cannot be generalized for other topologies. Indeed, it is essential to have 
the factor \i? x — y 2 and, therefore, Eq. (|107|) can be applied only once. However, we have now a 
result that can be used for the sunset diagram, without numerical instabilities, also at threshold. 
Numerical results will be shown in the next sections. 

For other topologies we will use the minimal BT-approach everywhere but at the corresponding 
normal thresholds and will adapt the present algorithm around them. Actually there are cases 
where the result is again very simple. 



- The S4 topology: 
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The S4 topology is given by 

7r 4 S , 4 = /i 2e / d n qi d n q 2 



(q\ + m?) ( (gi - q 2 f + ml) (<?f + m 
We easily derive the following relation: 



<?2 + p) + ml 




92 +P 



(148) 



Figure 8: The two-loop diagram S4. Arrows indicate the momentum flow. 
S4 = (e — 1) J dx S33 {x 2 p 2 \ mi, m 2 , M 3 



Ml 



where 



e - 



l)Sl^(x 2 p 2 ;m 1 ,m 2 ,M x 



2 2 1 / 2 1 2 2\ 1 2 

— p x + [p + m 4 — m 3 jx + m 3 , 



(149) 



/ dz 

Jo 



ln^(y,2;) + 



z- 1 



7 1 



C(2), 



+ InM 2 (2-lnM x 2 
(150) 



and ^ is given in Eq. ( |144j) , with p 2 — > x 2 p 2 . For other cases the result is more involved. 
- Other topologies: 

The key observation is that any two-loop diagram can be written as a mult i- dimensional integral, 
over external Feynman parameters x ext , of a kernel which is nothing but a generalized sunset 
topology with masses that depend on x ext , 



D 



n 

=1,N 



dx l ext (a, m 1 I 7, m 2 , p \ (3, m 3 ) . 



(151) 



Around thresholds, where nothing better will work, we make use of recurrence relations, for 
instance 

1 d 

(a + 1, mi I 7, m 2 , p \ (3, m 3 ) = — ^ (a, ra a 7, m 2 ,p \ (3, m 3 ) , (152) 

a am{ 

until we reach (1, mi | l,m 2 ,p \ 2,777,3) where a stable numerical result can be produced. At this 
point, despite the poor reputation enjoyed by numerical differentiation, we will perform all differ- 
entiations with respect to masses numerically. Only at the end the final integration over x ext will 
be executed[|. 

1 However, a detailed study of two-loop self-energies [^7| shows that, in these cases, numerical differentiation can 
be avoided in favor of alternative realizations of the raising procedure 
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3.1.6 Special cases and numerical results 



The whole algorithm presented in the previous sections has been translated into a FORM code pO 



and the output has been used to create a FORTRAN codeQthat computes S3. Some of the results 
are shown in the following sub-section where they are compared with the analytical results that 
are available in the literature for some special cases. When some of the masses are zero the result 
for S3 simplifies. For mi = we obtain 



£3 (p 2 ; 0, m 2 , m 3 J = m\ S 32 (p 2 ; 0, m 3 , m 2 ) + m\ S 33 (p 2 ; 0, m 2 , m 3 ) + S 3p (p 2 ; 0, m 3 , m 2 
For m 2 = or m 3 = we use Eq. ( |153| ) and the following properties: 

5*3 (p 2 ; m l7 m 2 , m 3 J = S 3 (p 2 ; m 2 , mi, m 3 J = 5*3 (p 2 ; m 3 , m 2 , m x 
Finally, if two masses are zero we use the appropriate permutation and 

S 3 (p 2 ; 0, m 2 , 0) = m\ S 32 (p 2 ; 0, 0, ra 2 J + S 3p (p 2 ; 0, 0, m 2 
For mi = m 3 = we have an analytical result |6l| (see also ref. []62]]). With 

S ■ c 

z = — + id, 
m A 



A uv - lnz, 



(153) 



(154) 



(155) 



(156) 



we derive the real and imaginary parts of S3 



ReS 3 (p 2 ;0,m,0) 



r2 1 



m 



1 



- + - ( 3 - -z - 2A UV + 2hiz) + - ( : - - ) In | : - I 

13 



1 



+ Re Li 2 (z) + 3 + - z A uv (z) - — z + A 2 uv (z) - 3 A uv (z) + - C(2) 



Im^^ 2 ; 0, m, 0) = m 2 Im 



z — 



In (!-*)+ Li, (2) 



(157) 



For the numerical evaluation we have used the NAG [33] Fortran library subroutine D01EAF that 
integrates a vector of similar functions (e.g. real and imaginary parts of the scalar diagram) over 
the same multi-dimensional hyper-rectangular region. In certain regions, e.g. around thresholds, 
one obtains moderately accurate results and for these cases we switch to the routine D01GDF, a 
multi-dimensional quadrature for general product regions, using number-theoretic methods. 

A detailed comparison is shown between the numerical approach based on the minimal BT- 
algorithm and the analytical result of Eq.( |157[ ) in Tab. |3[ There we use mi = m 3 = 0, m 2 = 
1.765 GeV and moreover we select e = 1 and u. = 1 GeV. Note that there is no loss of numerical 
precision near, but not too near, the normal thresholds. The requested precision for the integration 
routine is of a relative error of 10~ 5 . Comparison with the analytical result shows that the returned 
numerical error is overestimated with a very good agreement over a wide range of values for y^i. 
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16.333332 


-0.000196 




6572.7846 


-3862.1002 


2 


16.491314(184) 


-0.028545(7) 


100 


33332.656(159) 


-15629.648(102) 




16.491316 


-0.028545 




33332.642 


-15629.648 



Table 3: The sunset topology £3 for m\ = = and ni2 = 1.756 GeV as a function of 1/s. First entry 
is the minimal BT-approach, second entry is the analytical result of Eq.( |157| ). The UV-pole is e = 1 and 
the unit of mass is /i = 1 GeV. 



To investigate in more detail the situation around s = m\ we have produced Tab. [|. From it we 
see that by approaching y/s = rri2 up to X — ^fsjm^ = 3.4 x 10~ 3 we register an increasing numerical 
error on the real part, although the relative deviation from the analytical result is better than 10 -4 
signalling, once more, the conservative estimate of the numerical error. For the imaginary part, 
immediately after threshold, the numerical convergence is poor in a situation where the result is 
very small. At 1.4 x 10~ 2 away form threshold, however, we register again an excellent agreement. 

As explained in Section p.l.5| , however, around the normal threshold we can produce an accu- 
rate result by using Eq.( 147 ). Results for a scan around and at threshold are shown in Tab. |5|. 
Another comparison concerns the imaginary part of S3 which is finite in four dimensions and can be 
derived analytically by reduction to Legendre form of elliptic integrals. The complete expression is 
given in ref. |nj and, here, we only give a detailed comparison. The result is shown in Tab. |6| where 
we compare two different methods of numerical integration. We find better and better agreement 



2 FORTRAN may be archaic but never leaves you in solitude. 
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Table 4: A scan of the sunset, S3, topology for mi = 7713 = and m-2 = 1.756 GeV around the normal 
threshold s = (mi + m-2 + m,3) 2 . First entry is the minimal BT-approach, second entry is the analytical 
result of Eq. (|i~57 ). A * indicates poor numerical convergence in the numerical evaluation. Note, however, 



that this table is only for illustration and, in this region, one should always use Eq.( |147| ) (see Tab. |5|). 
The UV-pole is e = 1 and the unit of mass is (jl = 1 GeV. 

away from thresholds with a numerical error that, in general, overestimate the difference with the 
analytical result. Indeed, near threshold, we have a relative difference of 0.017%(0.005%) with a 
returned numerical error of 0.79%(0.74%). 

Finally, following the suggestion of ref. ||61| , we introduce a special combination: 



S c = S 3 (p 2 ; mi, m a , m 3 J - S 3 (p 2 ; m u 0, m 3 ) 

- S 3 (V; 0, m 2 , m 3 ) - S 3 (p 2 ; 0, 0, m 3 ) , (158) 



where infinite parts cancel. For small | p 2 |, essentially for | p 2 \< m 2 , the combination S c is known 
in the form of a multiple series in the variables 

Z\ — — z 2 — — 2, Z3 — — (15yJ 

VTi3 f«3 VTi3 

with s = — p 2 . The comparison is shown in Tab. |7] for p 2 space-like or time-like, showing excellent 
agreement with the analytical result derived from the expansion of a Lauriceila function. 

3.2 Tensor integrals 

Tensor integrals of the type 

(fii . . .fii I vi . . .Vj I l,m 1 I l,m2,2? I l,m 3 ) , (160) 
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Table 5: The same as in Tab. I but using Eq.( |147| ). Integration is performed with the NAG subroutine 
D01GDF. 

are usually decomposed into scalar integrals. For a general discussion we refer to the work of 
ref. |59| where the following functions are introduced: 



(e-l)Vl{ 2 = (p---P I P_-_P I 1, "7,1 I 1,^2, P I 2,777,3) = (i,P I j,P I 1, I l,m 2 ,p | 2,777,3), (161) 
« i 

where 2 + j < 3. It is easily shown that these integrals admit the following representation: 



112 



z(l - xj\ ^ y(l- yy /2 Q%{x, y, e) U n ~ e (x, y), (162) 



where U is given in Eq.( |102| ). Therefore, to each of these functions we apply the minimal BT- 
algorithm and derive expressions similar to those obtained for S33 = V^- The polynomials Q are 
as follows: 

Q™=2(ml + y y), Q° ° = - e -l, QT = 0, 



Q™ 1 = -2xyp\rni + yY), Qo =xyp 2 (l--) 1 Q 1U = 0, 
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Table 6: Im53, the imaginary part of the sunset topology, for m\ = m.3 = 80.448 GeV and m<i = 
91.1875 GeV as a function of \fs. First entry is the minimal BT-approach, first row integrated with 



D01EAF, second row with D01GDF and Eq.( 147 ); second entry is the analytical result of ref. |6l[| . 
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Table 7: ReS c (see Eq.fljggD for the definition) for m x = 10 GeV, m 2 = 20 GeV and m 3 = 100 GeV as a 
function of p 2 . First entry is the minimal BT- approach, second entry is the analytical result of ref. [61]. 
The UV-pole is e = 1 and the unit of mass is y, = 1 GeV. 
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(163) 



In our approach we avoid using derivatives with respect to the external parameters since the poor 
reputation enjoyed by numerical differentiation leads to search for techniques which avoid the 
explicit use of it. Conversely our numerical integration is usually performed with the help of the 
NAG JkJ Fortran library subroutine D01EAF that computes approximations to the integrals of a 
vector of similar functions over the same multi-dimensional hyper-rectangular region. Therefore 
all relevant integrals, connected to the 53 topology (any topology), are computed in one stroke 
without any noticeable loss of CPU time with respect to computing the scalar integral alone. 
Therefore we introduce additional tensor integrals of the following form: 
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where i,j = 0,1. They also can be cast into the form 

' r(e) r 1 , r 1 . r 



,pl,2|i 
' 112 
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i & J. 







dx / dy x(l-x) y(l-y) e / 2 g^(x,y,e)^(x,y), (165) 



with corresponding polynomials that will not presented in this paper. Other integrals can be easily 
derived by applying (a/57) partial integration, for instance 

1 



(hP I 3,P I 1 , m i I 1,^2, P I l,ms. 

v 7 e-l-(i+.?)/2 

? I j,p I 2, mi I l,m 2 ,p | l,m 3 ) + (p 2 + m 2 ,) | j,p | l,mi | 2,m 2 ,p | l,m 3 ) + 



m 
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™l (hP \j,p\ l, m i I ^, m 2,P | 2,m 3 ) + (i + l,p\ j,p \l,m 1 \ 2,m 2 ,p | l,m 3 ) + 

| j + l,p | l,mi | 2,m 2 ,p | l,m 3 ) . (166) 

The tensor integrals are easily computed and, once the method will be extended to cover all 
two-loop topologies, one may start organizing a realistic calculation. However, in any realistic 
calculation, the general strategy will be substantially different form the one usually adopted in 
one-loop exercises. It will be more convenient to write a Form code for any two-loop diagram G, 
or set of diagrams J2j Gj, that produces one expression of the form 

/ dxJ^Qi^V^ix), (167) 

and to apply the method directly to this integral instead of applying first the tensor decomposition 
and then to compute, separately, the various ij components of the total result. In a way one looses 
some attractive feature of the analytical one-loop approach, where few universal building blocks 
are needed and everything is algebraically reduced to those blocks. But, again, this is a typical 
example of a procedure motivated by our ability to derive analytical results. The perspective, 
here, is radically different, see ref. [|5(J for an epistemological discussion. 



4 Conclusions 



In this paper we have considered the Bernstein-Tkachov algorithm, proposed in |50|, which 



proves of some use in dealing with fast and accurate numerical evaluation of arbitrary multi-loop 
Feynman diagrams. 

As a prelude, we have applied the BT algorithm to one-loop diagrams with an arbitrary number 
of external legs, where reliable analytical results have been known for a long time, to show the 
feasibility of a project based on this new strategy and to discuss, with simple examples, the infrared 
divergent cases and the new approach that one will have to adopt in dealing with reduction of 
tensor integrals. 

After this brief excursus in the one-loop world we have considered the application of the 
algorithm to genuine multi-loop diagrams. Vacuum and tadpole diagrams and two-loop diagrams 
at zero external momentum allow for a full implementation of the method but, in general, we have 
not succeeded in deriving the basic ingredient, i.e. the polynomial P of Eq.([i])Q. Therefore, we 
have introduced a variant of the original proposal, a minimal BT-approach, which can be described 
as the convolution of the Ghinculov and Yao techniques, especially the use of complex Feynman 
parameters, with the basic BT-algorithm of raising powers in Feynman integrands. The main 
idea is related to the simple observation that we know how to apply the BT-iterative procedure 



3 To the best of our knowledge no example is known in the literature beyond the one-loop case of |5Q]. Working 
with some unrealistic case S. Uccirati has shown how to deal with cubic polynomials in two variables where, however, 



the coefficient B of Eq.(TJ) contains approximately 3 x 10 4 terms. 
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of raising powers for an arbitrary one-loop diagram. Therefore, given any two-loop diagram G 
we minimally apply BT functional relation to the one-loop sub-diagram L G G which has the 
largest number of internal lines. In this way the integrand can be made smooth, a part from the 
factor B of Eq.(|l|) which is now a polynomial in x s , the set of Feynman parameters needed for the 
complementary one-loop sub-diagram S G G with the smallest number of internal lines. Since 
the BT procedure does not introduce singularities through B, a part from the singularities of G 
itself in the external parameter space, before performing the x s -integration we move the contour 
into the complex hyper-plane, thus avoiding the crossing of apparent singularities. This is the 
procedure introduced in |59| . 

The minimal BT-approach has been applied to the sunset two-loop diagram and numerical 
results are shown, comparing them with the relevant literature. Whenever available, published 
results for the sunset show excellent agreement with our method which is fast and accurate in 
all regions, time-like or space-like external momentum, around normal thresholds and pseudo- 
thresholds. 

However, when we are close enough to the normal threshold the integration contour cannot be 
distorted and numerical instabilities appear. 

A preliminar analysis of other two-loop diagrams leads us to formulate the conjecture that these 
numerical instabilities are always connected with the leading Landau singularity of the original 
diagram. 

In these regions we have found more convenient to apply another algorithm for rasing powers 
which gives a very simple result, although with discontinuous imaginary part of the integrand. 
Therefore, very accurate predictions can be made also at the normal threshold. The extension of 
this variant of the rasing procedure to other two-loop topologies may require the introduction of 
numerical differentiation with respect to masses which is, in general, an unsatisfactory branch of 
numerical analysis. However, close to normal threshold, this alternative procedure represents the 
safest way for a numerical evaluation of diagrams, unless one want to adopt specific expansions 



that can be found in the literature. A preliminar analysis |57j has shown that, for two-point 
functions and one iteration, a transformation of the Feynman parameters can always be found 
that produces a coefficient B x s -independent. This B will vanish at some non- leading Landau 
singularity of the diagram where additional analytical work is needed before starting the numerical 
evaluation. 

Work is in progress to extend the method to all two-loop, two-(three-)point topologies and we 
plan to extend it to all two-loop diagrams with the evident goal of performing a complete two-loop 
analysis of the standard model predictions, in particular of quantities like sin 2 9 l eS at the level of 
1 x 10~ 6 theoretical precision. 
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5 Appendix A 



In this Appendix we give details about the analytical derivation of double and single UV poles 
for S33, defined in Eq . (|89|) . First we use the Laurent expansion of T(z) around z = 0, 



x — 



Next we introduce 



7 3 + 3C(2) + 2C(3) z 2 + 0(z 3 . 



x (1 — x)' 



and expand in e the various terms as follows: 
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Some of the integrals can be performed directly by using the following relations: 
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In this way we arrive at the following expression for S33: 
2 2 



C(2) 
n + 1' 



$33 = -^-- I 1 + 7 + In 7T + In — 1 + (7 + In vr) z 

+ ln4r (2 + 2 7 + In 4 ] + 2 lnvr (l+ln^) +27 
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-xy 3 (l -x) +yM 2 



(168) 
(169) 



(170) 



(171) 
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+ y dx j dyC 1 [2 + In (1 - y) + In (x - x 2 ) -2 lnf] [- ii/ 3 (l - 1) + i/M, 2 

— 4 / dx / dy y In £ ( 7 — In 7r — In -^r ] 
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4 2 

here we have introduced 



(172) 



M 2 = (1 — x) + x nl = x (1 — x) /i 2 ,. (173) 
next we restore £ = 2 (1 — x) x an d niake use of the relation Eg. (]107|) . After integration by parts, 



i.e. 



.1 yU Q 

Jo x oy 
the final result, Eqs.(|l06|)- ([108|) , follows. 



rl y n Q rl 

/ — — x = hi/i 2 -n / dyy n ' x \nx, 
Jo X oy Jo 



(174) 



6 Appendix B 

In this appendix we consider the case p 2 = —t, with t < 0. Since in this case 

U(x,y) = -y (1-y) t + ym\ + (1 - y) m 2 , (175) 
the quadratic form x-, f° r ^-channel diagrams is redefined as 

u = -t x , x = y (i-y) + y»l+ ( l ~v) /4 ( 176 ) 

with /i 2 = —vp? jt > 0. Furthermore, when raising the power /x in x^ f° r the t-channel case, we 
have to change coefficients accordingly: 



which gives 
and one should replace 
in the pole terms. 
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(177) 

(178) 
(179) 
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7 Appendix C 



In this Appendix we give the explicit integrand for S33, the other S^i terms following through 
the permutations of Eq . (|97|) . First we introduce auxiliary variables 



180) 



X = x(l-x), A = x(l-f4)+B, 
B = Xfi 2 x} £(x,y)=Xx(x,y)- 
The '+ '-distribution is defined, as usual, by its action on a generic test function g(x) 



dx 



g( x ) -0(1) f(x). 



dxg(x)f+(x) = 

JO 

We also use the notations 

ln+£(x,y) = \n£(x,y) - \n£(x, 1), L x = \nX, L y = \n(l-y). 
Then we derive the following expressions: 



S: 
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Double and single pole parts are defined in Eq. (|106|) and in Eq. (|108|) respectively. For the UV finite 
part we split the final result according to the dimension of the integrals over Feynman parameters, 



h — f dx I dyl 2 2 + [ dxL 

Jo Jo Jo 
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The result is as follows: 
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45 



V) {-<oABX + UAX 2 - 3A 2 X + A 2 B + 8BX 2 - 16X 2 ) ( ^p^p ) + 

+£(x, y) X {-4AB + AX(U - 2y) - A 2 (3 + 5y) + 245X(y + 1) - 16(y + 1)X 2 
-£(x, y) [2A + 16X(1 + y)] } ln + y) 



J 21 = -^ABX 2 + ±f^X 3 + UA 2 BX - ^A 2 X 2 



9 



13 



+X 



158 x4 

15 

586 



+AA 3 X X -4A 3 X B + y Bx3 

A (_ 20 £ 2 + ^BX + ^X 2 )+A 2 (9B - 72X) 

176 

y 



3 5 

+A 3 {l4X - AB^+X (24B 2 - UBX - — X 2 ) 



+.4X (^-25 2 + 2£X - AB + 5 AB 2 J In £(x, 0) + AfT In 0)L 
+ |^ (22XB 2 - 16X 2 £ - 126X 3 )+A 2 (-12X5 + 77X 2 - 5£ 2 ) 
+A 3 , (-15X + 5B^)-24B 2 X 2 + 405X 3 + 64X 4 ] ln£(x, 1) 



+ 



6AB 2 X - 30AX 3 - AA 2 BX + 17A 2 X 2 - A 2 B 2 

]n£(x,l)L x - A 2 B 2 ln 2 £(x,0) 



-3A 3 X X + A 3 X B 



W 2 X 2 + 8BX 3 + 16X 4 



+AX 



-6AB Z X + 30AX 6 + 4A Z BX - 17 A 2 X 2 + A A B 



|2 D 2 



+3A 3 X - AlB + 8£ 2 X 2 - 8£X 3 - 16X 4 



ln 2 £(x,l) 



In = y) (-24 + 16Xy) L y ln + £(x, y) + £(x, y) (lA - 40Xy) L x 

V) (-2A + 16Xy) ln£(x, y)L x + £(x, y) (-74 + 40Xy) ln£(x, y) 

y) (24 - mXy) In 2 y) + {-AX + AB - 2BX + 2X 2 ) ( ^^f ) 
-ln|(a;,y)> 



+£(*, y) (-24 + 4X) ( ln ^) ) y) [-24 + 4X (l + y)] ln + y) 



In = AX (-—AX + 24 2 -BX + —X 2 ) 
V 9 4 / 



11 



/ 133 ,n 110 o\ 

+4X ( AX - AAB + 64 2 + 9£X + X 2 L x 

V 6 fi / 



115 



+AX {—AX + AAB - 7 A 2 - 7BX - —X 2 ) ln£(x, 1 



133 



46 



+A (7 AX + AB- 2A 2 - 2BX - 6X 2 ) hx£(x, l)L x 
+AX (-1 AX - AB + 2A 2 + 2BX + 6X 2 ) In 2 f (ar, 1) 



'01 



13 1 

L x +]n£(x,l) I 00 = - + -((2) 



Similarly we derive the following expression for Sa p of Eq.(p9|): 
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where the single-pole part is given in Eq. (|142p . For the finite part we split again the 
contributions according to 
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The result for the Ji functions is as follows: 

J22 = £(x, y) \AXB{22Ay - 32) + A8A 2 X + 66A 2 B 
-64A 3 y - 200XB 2 - 128X 2 y#] L y ln+ £(x, y) 

+£ 2 (x, y) \AX{A2 - 192y) - 64^ 2 + 340X5 + mX 2 y\ L y ln+ £(x, y) 
-U0f(x,y)XL y ln + C(x,y) 

+ \AX(-524yB + 32B) - 80A 2 Xy - U8A 2 B + 139,4 3 y 
+528XB 2 + 256X 2 yB] L x £(x, y) 

+£ 2 (x, y) Ux(-56 + 484y) + 139A 2 - 924XS - 224X 2 y] L x + 396f(x, y)X L a 
+£(x, y) \AX{22AyB - 32B) + 48A 2 Xy + 66A 2 B - Q4A 3 y 
-200XB 2 - 128X 2 yB] ln£(x,y)L x 

+£ 2 (x, y) \AX{A2 - 192y) - 64A 2 + 340XS + 96X 2 y] Inf (sc, y)L x 
-140£ 3 (x,y)X ln£(x,y)L a 

+£(2, y) Ux(524yS - 32B) + 8QA 2 Xy + IA8A 2 B 
-139,4 3 y - 528XB 2 - 256X 2 yB] ln£(x,y) 



187) 
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+£ 2 (x, y) [AX(56 - 48%) - 139A 2 + 924X5 + 224X 2 y] ln£(x, y) 
-396£ 3 (x,y)X ln£(x,y) 

+£(x, y) \AX(-22AyB + 32£) - A%A 2 Xy - mA 2 B 
+UA 3 y + 200X5 2 + 128X 2 y#] In 2 £(x, y) 

+f(x, y) [AX"(-42 + 192y) + 64A 2 - 340XB - 96X 2 y] In 2 £(x, y) 
+140£ 3 (x,y)X ln 2 £(x,y) 



+£(z,y)[- 



4,4X5 + 4AX - A 2 X + A 2 B + 4X 2 5 - 4X" 



/ lng(x,y) \ 
V «-1 /- 



y-i 



+£(x, y) \AX{2AyB - AB) + 4AX 2 - A 2 X(1 + y) + 8A 2 5 
-6A 3 y - 32X£ 2 + 8X 2 £(1 + y) - 4X 3 (1 + y)] ln+ £(x, y) 
+£ 2 (x, y) [-22AXy - 6^ 2 + 52X5 - 4X 2 (1 + y)] ln+ £(x, y) 
-20£ 3 (x,y)X ln+£(x,y) 



J21 



59 



182 



31 



282 



— AX A B + ^AX i + — y^XB - — A 2 X 2 + ^X 4 B 



283 



75 



300 



75 



1234 
735" 



-X 4 



+ 



_20S 2 AX - -BAX 2 + — ^X 3 + -5A 2 X 
3 15 3 



2 v 2 



+ 



-— A 2 X 2 + 9B 2 A 2 + 5^ 3 X - 5BAl + 20-B^X 
5 

-4 In £(x, 0)A 2 .B 2 - In £(x, 0)L X A 2 B 2 
20AXB 2 + 20AX 2 J B - 40AX 3 - 20A 2 X£ 



412 
~2T 



-X 4 



_Zv 7 



+25,4 2 X 2 - 5A Z B Z - hA 6 X + hA 6 B - 20X*B* + 20X 



|2 D 2 



2 D 2 



+ 



ln£(x,l) 



4AXB 2 + 4AX 2 5 - 8AX 3 - 4A 2 XB 
+5A 2 X 2 - A 2 B 2 - ,4 3 X + A 3 B - AX 2 B 2 + 4X 4 j ln£(x, l)L a 
+£(x, O)^ 2 ^ 2 In 2 + [-4.4XB 2 - 4,4X 2 £ + 8.4X 3 + AA 2 XB 
-5A 2 X 2 + 4 2 B 2 + A 3 X - A 3 B + AX 2 B 2 - 4X 4 1 In 2 £(x, 1) 



J12 = £(x, y) [A(2 + 24y) - 16Xy - 30^] L y ln+ £(x, y) + 30£ 2 (x, y) L y ln+ £(x, y) 
+£(x, y) [A(-7 - 64y) + 40Xy + 86B] L x - 86£ 2 (x, y) 
+£(x, y) [A(2 + 24y) - 16Xy - 30#] ln£(x, y)L x + 30 ln£(x, y)L^ 2 (x, y) 
+£(x, y) \a(7 + 64y) - 40Xy - 86b] ln£(x, y) + 86£ 2 (x, y) ln£(x, y) 
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y) A(-2 - 24y) + 16Xy + 30B In 2 £(x, y) - 30£ 2 (x, y) In 2 £(x, y) 
'ka.£(x,y)' 



-A + 2X 



y-i 



y) + 3y) + 2X(1 + y)-6B ln+ £(x, y) + 6£ 2 (x, y) ln+ £(x, y) 



36 6 9 300 



+ 

+ 
+L 
+ 



12 2 3 15 

49 7 9 22 42 9 

—AX + -AB - 4A 2 XB X 2 

4 2 3 5 



3AX + AB-A 2 - 2XB - 2X 2 
3AX - AB + A 2 + 2XB + 2X 2 



ln£(x,l) 
ln 2 £(x,l) 



Joi = - [hLt(x, 1) - L X J J 00 = -— (188) 

8 Appendix D 

In this Appendix we derive results for the sunset diagram S 3 at p 2 = 0. In this case we can 
use the notation of rcf. 1631, i.e. 



(Mil, • • • , M ln | M 2i , . . . , M 2m | M 31 , . . . , M 3l 

1 



d n pd n q 17 ft TT 7 w 

l\ M fe =i P 2 + M 2 + M |. / (p _ g)2 + M | 



In these notations we have 



S 3 (0;m 1 ,m 2 ,m^= (mi | m 2 | ra 3 ). 



189) 



(190) 



Furthermore one easily derive 



(mi | 777,2 I 7713^= — (7 



-2 



;// I ///•> I //>:■()= — I ;//_> | ///;-, ) < " + (mi | WI2 | m^j ^ e 1 + (mi | 7772 | m 3) Q i (191) 



(-2) 



with singular parts given by 

( 
( 



mi 


m 2 


777,3) 


mi 


m 2 


m 3 ) 



(-2) 



2E 

1=1 



m„- 



(-1) 



mi I wi2 I m 3 



(-2) 



5^ m? [l -27-2 ln(vrm 2 ) 



(192) 



i=i 
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and a finite part that becomes 

(mi | mi | 7713 
3 



-( 



m\ m,2 mj, 



E 

i=i 



1 , , 

7 + in 7rm 

2 12 



)(-2) + ( mi I m2 I ms 
(7 + In 



7rm. — 1 



(-1) 



(193) 



Furthermore, the finite part is given in terms of the function 

Fi= dx Li 2 (l - fifj - - — — 2 In// 



ijX + bi(l — x) 
x(l — x^j 



^ 6- 

mf 



rat 



rn 



2 ■ 



(194) 



and 



3. 


k 


= 2 


for 


i = l, 


1- 


k 


= 3 


for 


i = 2, 


1. 


k 


= 2 


for 


i = 3. 



(195) 

Since p 2 = we can apply directly the original proposal of F. V. Tkachov | 5"0|] . First we combine 
the gi-dependent propagators with a Feynman parameter x and obtain 



<r gi 



(ql+mf) ((gi-g 2 ) + 



• 2-e/2 -p 
Z7T ' 1 



fife 







x(l — x) q\ + (1 — x) + x m 



-e/2 



(196) 



We can raise one power in the integrand and obtain 

,2 



V~ e / 2 = 4 



9 



[q 2 + ml — mf) +4 g 2 mf 



x 



{1 



2-e 



x — 



1 q 2 — m + m_ d 
q 2 dx 



•e/2 



where ra± = mi i m 2 and 

V = x(l — x)q + (1 — x) ra 1 + xm 2 = x(l — x) ( q + m 



Next we observe that 



(q 2 + m 2 2 - '+Aq 2 m\ = (q 2 + ra 2 ^ (q 2 + 
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m 



(197) 



'198) 



(199) 



and perform the remaining g-integration which requires the introduction of two additional Feyn- 
man parameters, y and z. This procedure will introduce two monomials in y, z, 



U±(y, z) = (m| - mfj y + (m\ - m|) z + m 2 x , 
which appear with power — 1 — e, so that we may use the following identity: 



mi 



l+ y -dy + -d Z 

e e 



(200) 



(201) 



After that we proceed as usual, namely we integrate by parts and expand in Laurent series around 
e = 0. The final result is 



2 2 

mi — m_ 



S , f n (0;mi,m 2 ,m 3 ) = 

+ / dx dy 
Jo Jo 



f dx f dy I dz K 3 
Jo Jo Jo 



K 2 + / dxK x +K . 
Jo J 



(202) 



We now introduce auxiliary variables 

W±(x, y, z) = x(l — x^j zm\ + (y — z}m\ + - (l — yj (m + + m^j —Am + m_x 

X± = xm\ + (l — x^j m 2 ± , A = (m + + m_) —4m + m_x, X = x (l — x^j, 



L y = \n(l-y), L x = \nX, 



(203) 



and derive 



K 



. .; = m\ K 2 ^ + - m 2 _ K 2 , 3 -, 



K 3± = \nW ± (x,y,z) {3(L X -T ±Jy) 



A - l2W±{x, y, z)\ ~2 A ~ 54W±(x, y, z) 
3 ,\nW±(x,y,zy 



-3 In W ± (x, y, z) [A - 12W ± (x, y, z)] } - °-A (^^iiZ) 
K 2± = l -A In W±(x, 1, y){3 L y + 1 + 3 [L x - In W ± (x, l,y)]} 



3 
1' 
3 



K\ = —m 2 (m + m^ H — m\ H — m 2 _) In — lnm? 

1 4 " X+ 



+m 2 (m+rri- m 2 , m 2 _ ) In — lnm? 



1 



m + m_ ln% (m+ — m 2 _) \nm\ — lnm 2 . 
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+ 



3m + m_ + (m 2 , + m 2 _) 



(m + — m_) In x+ In — 2 

m 2 



1 r 
4 

- mi) lnx+(l - - lnx+) - ^m 2 _(3m 2 + + m 2 __ 
X- 



x(ln^_-ln'x + -21n^- 



X 4 



K = +- 



+ 



ml + i« + ml)] ((2) - i 
m + m_ — 2~( m + + m2 ) 



m + m_ + — (m + + m. 



lnm^ 



lnmj + 



4 



in m 2 + —m^ r( m 4- + 



17. 

T 



(204) 



The above results have been coded in a FORTRAN program and compared with those obtained 
with the analytical result of [64|. The comparison is shown in Tab. || where we use e = 1 and also 
the jj, — 1 GeV, where /x is the unit of mass. 



mi [GeV] 


m 2 [GeV] 


m 3 [GeV] 


BT-numerical 


ref. 


54 




3 


2 


4 


147.63632 


147.63634 


6 


2 


4 


491.36054 


491.36060 


9 


2 


4 


1437.2692 


1437.2693 


90 


2 


4 


585240.60 


585240.63 



Table 8: The sunset topology S3 for p 2 = 0. First entry is the numerical BT-approach, second entry 
is the analytical result of Eqs.( |193| )-( |194] ) (see ref. |54|). The UV-pole is e = 1 and the unit of mass is 
fi = 1 GeV. 



52 



References 

[1] A. Pais, Oxford, Uk: Clarendon ( 1986) 666 P. New York, Usa: Oxford Univ. Pr. ( 1986) 
666p. 

[2] P. Mastrolia and E. Remiddi, Nucl. Phys. Proc. Suppl. 89 (2000) 76. 

[3] K. G. Chetyrkin and A. Retey, Nucl. Phys. B 583 (2000) 3 |[hep-ph/9910332|1 . 

[4] K. G. Chetyrkin and A. Retey, |hep-ph/0007088| 

[5] G. 't Hooft and M. Veltman, Nucl. Phys. B 153 (1979) 365. 



[6] Z. Bern, L. Dixon and D. A. Kosower, Nucl. Phys. B 412 (1994) 751 [[hep-ph/9306240 



[7] Z. Bern, L. Dixon, D. C. Dunbar and D. A. Kosower, Nucl. Phys. B 425 (1994) 217 flhep- 
ph/940322^. 



[8] Z. Bern, L. Dixon, D. C. Dunbar and D. A. Kosower, Nucl. Phys. Proc. Suppl. 39BC (1995) 
146 Khep-ph/94092141 ; 



Z. Bern, L. Dixon and D. A. Kosower, Ann. Rev. Nucl. Part. Sci. 46 (1996) 109 ||hcp 



ph/9602280 



Z. Bern, L. Dixon, D. C. Dunbar and D. A. Kosower, [hep-ph/9706447 



[9] G. Passarino and M. Veltman, Nucl. Phys. B 160 (1979) 151. 
[10] J. Fleischer, F. Jegerlehner and O. V. Tarasov, Nucl. Phys. B 566 (2000) 423 jxep- 



ph/9907327 



[11] O. V. Tarasov, Acta Phys. Polon. B 29 (1998) 2655 ||hep-ph/9812"250| 
[12] M. J. Veltman and D. N. Williams, |hep-ph/9306228[ 



[13] A. A. Vladimirov and D. V. Shirkov, Sov. Phys. Usp. 22 (1979) 860 [Usp. Fiz. Nauk 129 
(1979) 407]. 

[14] S. G. Gorishnii, S. A. Larin, L. R. Surguladze and F. V. Tkachov, Comput. Phys. Commun. 
55 (1989) 381. 

[15] F. V. Tkachov, New methods for evaluation of multi-loop Feynman diagrams, PhD thesis, 
INR, Moscow, March 1984. 

[16] F. V. Tkachov, Phys. Lett. B 100 (1981) 65. 

[17] K. G. Chetyrkin and F. V. Tkachov, Nucl. Phys. B 192 (1981) 159. 

[18] F. V. Tkachov, Theor. Math. Phys. 56 (1983) 866 [Teor. Mat. Fiz. 56 (1983) 350]. 

53 



[19] D. J. Broadhurst and A. G. Grozin, [hep-ph/9504400 



[20] V. V. Belokurov and N. I. Usyukina, Theor. Math. Phys. 41 (1979) 945 [Teor. Mat. Fiz. 41 
(1979) 157]. 

[21] D. I. Kazakov and A. V. Kotikov, Theor. Math. Phys. 73, 1264 (1988) [Teor. Mat. Fiz. 73, 
348 (1988)]. 

[22] N. I. Usyukina and A. I. Davydychev, Phys. Lett. B 298, 363 (1993). 
[23] V. A. Smirnov, Phys. Lett. B 500 (2001) 330 [|hep-ph/0011056| . 
[24] A. V. Kotikov, Phys. Lett. B 254 (1991) 158. 



[25] E. Remiddi, Nuovo Cim. A 110 (1997) 1435 |hep-th/9711188 



M. Caffo, H. Czyz, S. Laporta and E. Remiddi, Nuovo Cim. A 111 (1998) 365 |pej> 
th/9805118|1 ; 



M. Caffo, H. Czyz, S. Laporta and E. Remiddi, Acta Phys. Polon. B 29 (1998) 2627 JEep- 
th/9807119|| ; 



T. Gehrmann and E. Remiddi, Nucl. Phys. B 580 (2000) 485 |hep-ph/99l2"H2^| 1; 
M. Caffo, H. Czyz and E. Remiddi, Nucl. Phys. B 581 (2000) 274 [|hep-ph/9912501 



[26 
[27 
[28 
[29 
[30 
[31 
[32 
[33 
[34; 



T. Gehrmann and E. Remiddi, Nucl. Phys. Proc. Suppl. 89 (2000) 251 ||hep-ph/ 0005232 
T. Gehrmann and E. Remiddi, [hep-ph /0008287| ; 
T. Gehrmann and E. Remiddi, [hep-ph/0101147| ; 
M. Caffo, H. Czyz and E. Remiddi, |hep-ph/0103Qll . 

S. Laporta, Int. J. Mod. Phys. A 15 (2000) 5087 ||hep-ph/0 102033] . 

O. V. Tarasov, Nucl. Phys. Proc. Suppl. 89 (2000) 237. 

O. V. Tarasov, Nucl. Phys. B 502 (1997) 455 ||hep-ph/97033T9| . 

F. V. Tkachov, Phys. Lett. B 124 (1983) 212. 

S. G. Gorishnii, S. A. Larin and F. V. Tkachov, Phys. Lett. B 124 (1983) 217. 

G. B. Pivovarov and F. V. Tkachov, In *Tbilisi 1986, Proceedings, Quarks '86* 227-234. ■ 
F. V. Tkachov, Phys. Atom. Nucl. 56 (1993) 1558 [Yad. Fiz. 56 (1993) 180]. 

F. V. Tkachov, Phys. Lett. B 412, 350 (1997) [ |hep-ph/9703424[ . 

J. Fleischer, V. A. Smirnov, A. Frink, J. G. Korner, D. Kreimer, K. Schilcher and J. B. Tausk, 
Eur. Phys. J. C 2 (1998) 747 ||hep-ph/9704333l ; 

M. Beneke and V. A. Smirnov, Nucl. Phys. B 522 (1998) 321 [|hep-ph/9711391|| ; 

V. A. Smirnov and E. R. Rakhmetov, Theor. Math. Phys. 120 (1999) 870 [Teor. Mat. Fiz. 



54 



120 (1999) 64] [|hep-ph/9812529|1 ; 

A. I. Davydychev and V. A. Smirnov, Nucl. Phys. B 554 (1999) 391 ||hep-ph/ 9903328 
V. A. Smirnov, Phys. Lett. B 460 (1999) 397 [|hep-ph/9905323|l ; 



V. A. Smirnov and O. L. Veretin, Nucl. Phys. B 566 (2000) 469 ||hep-ph/ 990738^ ; 
V. A. Smirnov, Phys. Lett. B 465 (1999) 226 |[hep-ph/990747ll 



P. A. Baikov and V. A. Smirnov, Phys. Lett. B 477 (2000) 367 [|hep-ph/0001192|| ; 
V. A. Smirnov, Phys. Lett. B 491 (2000) 130 | |hep-ph/0007032] ; 



[36] V. A. Smirnov, |hep-ph/010lT52 



[37] F. V. Tkachov and K. G. Chetyrkin, In * Zvenigorod 1979, Proceedings, Group Theoretical 
Methods In Physics, Vol. 2*, 223-227. 

[38] K. G. Chetyrkin, A. L. Kataev and F. V. Tkachov, Nucl. Phys. B 174 (1980) 345. 

G. B. Pivovarov and F. V. Tkachov, Int. J. Mod. Phys. A 8 (1993) 2241 ||hep-ph/9612287 



[40] J. Fleischer, O. V. Tarasov and M. Tentyukov, Nucl. Phys. Proc. Suppl. 89 (2000) 112. 
[41] O. V. Tarasov, Nucl. Phys. B 480 (1996) 397 ||hep-ph/9606238l . 



[42] O. V. Tarasov, Phys. Rev. D 54 (1996) 6479 ||hep-th/96060lS|l . 

[43] J. Fleischer and O. V. Tarasov, Nucl. Phys. Proc. Suppl. 37B (1994) 115 ||hep-ph/ 940 7235 



[44] D. J. Broadhurst, Phys. Lett. B 164 (1985) 356; 
D. J. Broadhurst, Z. Phys. C 47 (1990) 115; 
D. J. Broadhurst, Z. Phys. C 54 (1992) 599; 

D. J. Broadhurst, J. Fleischer and O. V. Tarasov, Z. Phys. C 60 (1993) 287 [|iep-ph/9304303 | 



D. J. Broadhurst and D. Kreimer, Phys. Lett. B 426 (1998) 339 [ |hep-th/96120TT 



[45] G. Degrassi, P. Gambino and A. Vicini, Phys. Lett. B 383 (1996) 219 [|iep-ph/9603374|| ; 
G. Degrassi, P. Gambino and A. Sirlin, Phys. Lett. B 394 (1997) 188 ||hep-ph/961i363| ; 
G. Degrassi and P. Gambino, Nucl. Phys. B 567 (2000) 3 ||hep-ph/ 9905472] . 

[46] A. Freitas, S. Heinemeyer, W. Hollik, W. Walter and G. Weiglein, Nucl. Phys. Proc. Suppl. 
89 (2000) 82 ||hep-ph/ 00071291 



A. Freitas, W. Hollik, W. Walter and G. Weiglein, Phys. Lett. B 495 (2000) 338 Fep- 
ph/0007091|| ; 



A. Freitas, S. Heinemeyer, W. Hollik, W. Walter and G. Weiglein, |hep-ph/0101260 
S. Heinemeyer and G. Weiglein, [hep-ph/0102317 . 



55 



[47] T. van Ritbergen and R. G. Stuart, Phys. Rev. Lett. 82 (1999) 488 [|hep-ph/9808283| . 

[48] F. A. Berends, A. I. Davydychev and N. I. Ussyukina, Phys. Lett. B 426 (1998) 95 |iep- 



[49 

[50 

[51 
[52 
[53 
[54 

[55 
[56 

[57; 

[58 
[59 



ph/9712209 



S. Groote, J. G. Korner and A. A. Pivovarov, Phys. Lett. B 443 (1998) 269 [|liep-ph/ 9805224 



S. Groote, J. G. Korner and A. A. Pivovarov, Nucl. Phys. B 542 (1999) 515 [|hep-ph/9806402 



A. I. Davydychev and V. A. Smirnov, Nucl. Phys. B 554 (1999) 391 ||hep-ph/9903"328^ ; 

S. Groote, J. G. Korner and A. A. Pivovarov, Eur. Phys. J. C 11 (1999) 279 [|iep-ph/9903412| ; 

A. Bashir, R. Delbourgo and M. L. Roberts, |hep-th/010TT48| 



M. Caffo, H . Czyz, S. Laporta and E. Remiddi, Nuovo Cim. A 111 (1998) 365 |pfij> 
th/9805118|| ; 



M. Caffo, H. Czyz and E. Remiddi, Nucl. Phys. B 581 (2000) 274 [|hep-ph/9912501 
M. Caffo, H. Czyz and E. Remiddi, |hep-ph/0103014 . 

R. J. Eden, P. V. Landshoff, D. I. Olive, and J. C. Polkinghorne, The analytic S matrix, 
Cambridge Univ. Press, 1966. 287p. 

F. V. Tkachov, Nucl. Instrum. Meth. A 389 (1997) 309 |[hep-ph/9609429fl ; 
L. N. Bertstein, Functional Analysis and its Applications, 6(1972)66. 

D. Y. Bardin, L. V. Kalinovskaya and F. V. Tkachov, |hep-ph/0012209| . 

G. Passarino, |hep-ph/0101"299 . 



D. Bardin and G. Passarino, Oxford, UK: Clarendon (1999) 685 p. 

F. V. Tkachov, Int. J. Mod. Phys. A 8 (1993) 2047 ||hep-ph/96l2"2S4 
F. V. Tkachov, Sov. J. Part. Nucl. 25 (1994) 649 [|hep-ph/9701272|| . 



G. Devaraj and R. G. Stuart, Nucl. Phys. B 519 (1998) 483 [|hep-ph/9704308 



P. Cvitanovic and T. Kinoshita, Phys. Rev. D 10 (1974) 3978. 

G. Passarino and S. Uccirati, Numerical Evaluation of Two-Loop Self-Energies, work in prepa- 
ration. 

G. 't Hooft and M. Veltman, Nucl. Phys. B 44 (1972) 189. 

A. Ghinculov and Y. Yao, Phys. Rev. D 63 (2001) 054510 ||hep-ph/00063T4] ; 
A. Ghinculov and Y. Yao, Mod. Phys. Lett. A 15 (2000) 1967 ||hep-ph/ 0004201 



A. Ghinculov and Y. Yao, Mod. Phys. Lett. A 15 (2000) 925 [|hcp-ph/0002211|] ; 
A. Ghinculov and Y. Yao, Nucl. Phys. B 516 (1998) 385 ||hep-ph/9702"26lj 



[60] J. A. Vermaseren, |math-ph/0010025| . 



56 



[61] F. A. Berends, A. I. Davydychev and V. A. Smirnov, Nucl. Phys. Proc. Suppl. 51C (1996) 
301 ||hep-ph/960624l ; 



F. A. Berends, A. I. Davydychev and V. A. Smirnov, Nucl. Phys. B 478 (1996) 59 [[hep- 
ph/9602396|| 



S. Bauberger, F. A. Berends, M. Bohm and M. Buza, Nucl. Phys. B 434 (1995) 383 [fiep^ 
ph/9409388| 



S. Bauberger, M. Bohm, G. Weiglein, F. A. Berends and M. Buza, Nucl. Phys. Proc. Suppl. 
37B (1994) 95 ||hep-ph/ 94064041 ; 



F. A. Berends and J. B. Tausk, Nucl. Phys. B 421 (1994) 456. 
[62] M. Caffo, H. Czyz, S. Laporta and E. Remiddi, Nuovo Cim. A 111 (1998) 365 flh"ei> 



th/9805118| 



[63] NAG Fortran Library, Mark 19, The Numerical Algorithms Group Ltd, Oxford UK. 1999. 
[64] J. van der Bij and M. Veltman, Nucl. Phys. B 231 (1984) 205. 



57 



Imx 




Figure 9: The complex a>plane with the cut of Inx (Eq.(|34|)) internal to the interval [aC, a;Z] (Eq. (|119| )), 
Dashed lines indicate the sign of the imaginary part of the logarithm when approaching the real axis on 
the cut. An example of the deformation of the integration path is shown. 
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Figure 10: The complex x-plane with the cut of lnx (Eq.(|34])) external to the interval [aC,a;+] (Eq. (|119| )). 
Dashed lines indicate the sign of the imaginary part of the logarithm when approaching the real axis on 
the cut. An example of the deformation of the integration path is shown. 
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Figure 11: The behavior of x L and x R , the roots of the quadratic of Eq.(144)) as a function of y for the 
sunset diagram at the normal threshold s = (mi + m% + m^) 2 . The figure shows that the integration 
contour is pinched by the two branch points x L)R which are the roots of Eq.(144). 
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Figure 12: The behavior of x L and x R , the roots of the quadratic of Eq. (|144j) ) as a function of 
sunset diagram at the pseudo threshold s = (m\ + m% — raj) 2 . 
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